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Exact  unconditional  hypothesis  testing  methods  in  the  context  of  the  2x2  bino- 
mial trial,  2x2  multinomial  trial,  and  the  comparison  of  correlated  proportions  ex- 
hibit a  superior  power  advantage  when  compared  with  exact  conditional  approaches. 
However,  until  now  there  have  been  no  attempts  made  to  compare  unconditional  and 
conditional  methods  in  a  more  elaborate  setting,  i.e.  testing  scenarios  that  involve 
more  than  just  a  single  nuisance  parameter.  This  was  due  to  the  unavailability  of 
techniques  present  for  computing  the  maximum  of  null  power  functions  involving  sev- 
eral nuisance  parameters. 

A  methodology  is  developed  herein  which  provides  a  means  to  compute  the  max- 
imum of  null  power  functions  that  are  of  two  possible  forms.  These  two  forms  cor- 
respond to  multinomial  and  binomial  likelihoods,  which  can  contain  several  nuisance 
parameters.  As  special  cases,  the  2x3  contingency  table  to  compare  two  independent 
trinomial  distributions,  and  the  2  x  2  x  2  contingency  table  to  compare  independent 
proportions  at  two  different  strata  are  considered. 

vii 


For  the  equal  sample  size  case,  exact  critical  values  of  the  Pearson  chi-square  test 
for  comparing  two  independent  trinomial  distributions  are  computed  and  tabulated 
for  n  =  10(1)70,  a  =  0.025  and  a  =  0.05.  To  facilitate  comparisons,  a  technique 
for  determining  sample  size  requirements  when  using  the  conditional  approach  was 
developed.  For  a  =  0.025  and  a  power  of  0.80,  sample  size  requirements  for  the 
conditional  test  were  often  times  smaller  than  those  for  the  exact  unconditional  chi 
square  and,  infrequently,  the  sample  size  requirements  were  larger.  Graphs  of  the 
null  power  functions  and  the  comparison  of  the  critical  regions  suggest  that  the  con- 
ditional p-value  is  a  better  test  statistic  than  Pearson's. 

For  the  2  x  2  x  2  table  exact  critical  values  of  the  Mantel  Haenszel  test  statistic 
with  continuity  correction  are  computed  and  tabulated  for  n  =  10(1)50,  a  =  0.025 
and  a  =  0.05,  where  the  sample  size  n  is  the  same  for  each  binomial  trial.  Again, 
a  method  was  developed  to  determine  sample  size  requirements  for  the  conditional 
approach,  so  comparisons  could  be  made.  For  a  —  0.025  and  a  power  of  0.80,  sample 
size  requirements  for  each  method  were  established.  The  results  proved  inconclusive, 
as  neither  method  outperformed  the  other.  However,  graphs  of  the  null  power  func- 
tions indicated  that  the  conditional  p-value  is  much  more  stable  than  the  Mantel 
Haenszel  statistic,  and  comparisons  of  the  respective  critical  regions  indicate  that  the 
conditional  p-value  is  the  better  statistic. 
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CHAPTER  1 
LITERATURE  REVIEW 

1.1  Introduction 

In  hypothesis  testing  procedures  where  both  the  null  and  alternative  hypotheses 
are  simple,  the  optimal  solution  for  a  given  level  a  is  completely  specified  by  the  Ney- 
man  Pearson  lemma.  However,  when  the  distribution  of  the  test  statistic,  under  H0, 
depends  on  additional  parameters  we  are  not  interested  in,  formally  called  nuisance 
parameters,  then  the  solution  is  anything  but  clear.  Suppose,  for  example,  that  we 
are  interested  in  conducting  a  hypothesis  test  about  the  parameter  9,  but  the  distri- 
bution of  the  statistic  used  for  the  testing  procedure  also  depends  on  the  nuisance 
parameter  7.  The  problem  we  are  faced  with  is  how  to  conduct  our  hypothesis  test 
in  the  presence  of  this  unknown  nuisance  parameter. 

There  are  numerous  solutions  to  the  nuisance  parameter  problem,  and  Basu  (1977) 
summarizes  these.  Of  the  methods  he  reviews,  only  three  will  be  considered  in  what 
follows.  It  also  should  be  noted  that  of  these  three,  only  two  are  regularly  used  in 
practice.  The  first,  and  perhaps  most  common  method,  is  using  asymptotic  theory. 
Here,  all  parameters  except  those  being  tested  are  replaced  by  estimates.  In  most 
cases  the  estimates  are  the  maximum  likelihood  estimates,  which  under  certain  reg- 
ularity conditions  are  consistent.  Since  they  are  consistent,  we  can  use  limit  theory 
results  such  as  Slutsky's  theorem  (in  conjunction  with  the  central  limit  theorem)  to 
derive  the  asymptotic  null  distribution  of  the  test  statistic.  It  frequently  turns  out 
that  the  resulting  distribution  is  simple  and  well  tabulated.  As  an  example,  the  Pear- 
son chi-square  test  statistic  for  comparing  two  independent  binomial  proportions  has 
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an  asymptotic  chi-square  null  distribution  with  one  degree  of  freedom.  The  advantage 
of  using  asymptotic  methods  is  simplicity.  However,  as  this  is  an  asymptotic  result, 
it  may  not  work  well  for  small  samples. 

The  second  method  is  the  conditional  approach.  Since  we  often  work  with  likeli- 
hoods that  are  members  of  exponential  families,  we  can  use  the  methods  of  Lehmann 
(1982)  to  eliminate  nuisance  parameters.  For  example,  suppose  we  are  interested  in 
testing  a  hypothesis  concerning  the  parameter  9  in  the  presence  of  an  unknown  nui- 
sance parameter  7.  Also,  suppose  that  X  and  Y  are  random  variables  that  are  jointly 
sufficient  for  (#,7),  and  the  joint  probability  distribution  of  (X,Y)  is  a  member  of  an 
exponential  family,  i.e. 

fxY{x,y;0,i)  =  C{6,  -y)exp(tx(8)x  +  t2(y)y  +  g(x,y)). 

Since  7  is  unknown,  we  can  restrict  our  class  of  tests  to  a-similar  tests.  Using  the 
conditional  distribution  of  X  given  Y,  we  can  develop  a  uniformly  most  powerful  un- 
biased (UMPU)  test  for  the  hypothesis  about  9.  The  major  criticism  of  this  method 
is  that  if  X  and  Y  are  discrete,  the  conditional  distribution  of  X  given  Y  will  be 
discrete  as  well,  so  that  a  randomization  procedure  is  needed  to  obtain  a  test  size 
of  exactly  a.  However,  since  randomized  tests  are  not  used  in  practice,  the  result- 
ing non-randomized  conditional  test  will  be  conservative.  Another  disadvantage  of 
the  conditional  method  is  that  we  are  restricted  to  a  conditional  sample  space.  All 
statistical  inference  procedures  will  be  done  with  respect  to  this  conditional  space. 
This  will  pose  a  problem  to  the  non-statisticians,  in  that  the  interpretation  of  the 
attained  significance  level  (p-value)  is  no  longer  appealing,  or  for  that  matter  eas- 
ily understandable.  When  conditional  methods  are  used  in  the  applied  literature, 
such  as  medicine  or  agriculture,  the  conditional  nature  of  the  p-value  is  almost  never 
mentioned  in  the  publication. 

The  third  technique  is  the  maximization  method.  Here,  the  size  of  a  test  is 
determined  by  maximizing  the  null  power  function  over  the  domain  of  the  nuisance 
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parameter  space.  In  this  way  we  are  able  to  eliminate  the  nuisance  parameters  and 
determine  the  form  of  the  critical  region  for  the  test  statistic  we  are  considering. 
For  example,  suppose  we  collect  n  data  points,  x  =  (x\,X2,  ...,£n)j  and  form  a  test 
statistic,  T(x),  which  supports  Ha  if  it  is  large,  i.e.  T(x)  >  a*.  Suppose  further 
that  when  H0  is  true  the  probability  distribution  of  T(x)  depends  on  the  unknown 
parameter  6.  For  a  test  of  size  a,  a*  would  be  determined  by  solving 

a*  =  inf{a  €  9?  :  sup7r(0)  <  a} 

where  n(9)  is  the  null  power  function,  which  depends  on  the  unknown  parameter  9. 
This  method  caters  to  the  worst  possible  configuration  of  the  nuisance  parameters, 
ensuring  that  the  type  I  error  rate  will  never  exceed  a.  An  advantage  of  this  method 
is  that  it  is  easily  explained  to  non-statisticians,  since  for  a  given  value  of  the  nui- 
sance parameter  they  are  familiar  with  the  methods.  Letting  the  nuisance  parameter 
vary  should  not  pose  any  problem,  and  they  should  be  able  to  develop  an  intuitive 
feeling  for  the  method.  Other  advantages  will  be  demonstrated  in  what  follows.  The 
drawback  to  this  method  is  the  complexity  of  maximizing  the  null  power  function. 

The  maximization  method  has  not  been  extensively  investigated  in  the  literature. 
However,  in  the  cases  where  it  was  studied,  the  results  were  favorable.  With  the 
exception  of  Mehta  and  Hilton  (1993)  and  Berger  and  Boos  (1994),  all  of  the  research 
pertaining  to  this  method  has  been  with  respect  to  the  2  x  2  contingency  table.  The 
reason  for  this  is  the  computational  complexity  involved  in  maximizing  the  null  power 
function  associated  with  tables  larger  than  2x2,  and  indeed  even  the  2  x  2  case  is 
quite  complicated.  However,  since  the  results  for  the  2x2  case  were  so  favorable 
(as  will  be  demonstrated),  it  would  seem  beneficial  to  develop  a  method  that  would 
handle  these  computational  difficulties  in  larger  contingency  tables.  We  develop  such 
a  method  in  chapter  2,  and  demonstrate  its  utility  in  two  applications.  The  first, 
presented  in  chapter  3,  is  the  2  x  2  table,  where  our  interest  is  the  comparison  of 
two  independent  trinomial  random  variables.  The  second  application  pertains  to  the 
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2x2x2  contingency  table  for  comparing  independent  proportions  across  two  strata, 
and  is  given  in  chapter  4. 

1.2    A  History  of  the  Maximization  Method 

The  maximization  method  was  first  suggested  by  Barnard  (1945),  with  details 
given  in  Barnard  (1947).  He  developed  a  method  for  the  analysis  of  the  2x2  com- 
parative trial,  which  provided  an  alternative  to  the  conditional  method  of  Fisher 
(i.e.  Fisher's  exact  test).  Although  he  mentioned  an  underlying  hypergeometric  and 
multinomial  structure  for  the  table,  the  main  focus  was  on  binomial  sampling.  In 
particular,  assume  we  have  two  independent  binomial  random  variables  X\  and  Xi 
with  parameters  (m,  6\)  and  (n,  62)  respectively.  It  is  of  interest  to  test  the  hypothesis 
H0  :  d\  —  62  versus  Ha  :  6\  ^  62.  For  small  sample  sizes,  the  usual  method  at  the 
time,  and  presently  as  well,  was  Fisher's  exact  test.  This  test  conditions  on  the  total 
number  of  successes,  and  sums  the  appropriate  hypergeometric  probabilities  in  order 
to  obtain  an  observed  significance  level  (p- value).  For  example,  suppose  there  are  t 
successes  observed.  Let 


H{xut)  = 


max(0,  t  —  n)  <  X\  <  min(£,  m). 


Then  the  attained  significance  level  is  given  by 

c 

where  C  =  {y  :  H(y,t)  <  H(xi,t)}.  Although  this  provides  a  very  simple  way 
of  eliminating  the  nuisance  parameter  8,  the  common  value  of  9X  and  62  under  H0, 
this  conditional  distribution  is  highly  discrete  for  small  sample  sizes.  The  result  is 
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a  conservative  testing  procedure,  which  is  what  motivated  Barnard  to  develop  the 
maximization  method. 

For  testing  HQ  :  0i  =  92  versus  Ha  :  9\  ^  92,  Barnard  believed  that  the  form  of 
the  rejection  region  (R)  should  possess  two  characteristics  independent  of  probability 
considerations.  These  were  called  symmetry  (S)  and  convexity  (C),  and  are  defined 
as  follows 

1.  Symmetry  (S)  :  if  (xi,x2)  G  R,  then  we  must  have  that  (m  -  Xi,  n  -  x2)  G  R. 
This  is  logical,  since  for  the  two-sided  case,  testing  H0  :  9\  =  92  is  equivalent  to 
testing  J?o  :  1  -  0i  =  1  -  02- 

2.  Convexity  (C)  :  if  (xi,x2)  G  R  and  if  {xux2)  lies  above  the  graph  of  x2  =  ^xt, 
then  any  point  "north"  or  "west"  of  (xi,x2)  should  be  a  member  of  R.  Similarly, 
if  (xx,x2)  G  R  and  if  (x\,x2)  lies  below  the  graph  of  x2  —         then  any  point 
"south"  or  "east"  of  (xi,x2)  should  be  a  member  of  R. 

To  determine  the  form  of  R,  the  above  conditions,  S  and  C,  are  used  in  conjunction 
with  the  following  condition 

3.  The  maximum  of  the  null  power  function  over  the  domain  of  the  nuisance 
parameter  space  should  not  exceed  a,  the  significance  level  of  the  test,  i.e. 


r  V  Xx  /  V  x'2  ) 

The  rejection  region  R  that  is  formed  on  the  basis  of  these  three  criteria  is  not 
necessarily  unique.  In  order  to  attain  the  largest  power,  the  region  which  has  the 


max  7t(0)  <  a 

O<0<1  — 


where  9  is  the  common  value  of  9i  and  92  under  H0,  and 


6 


greatest  number  of  points  is  selected.  The  fact  that  this  region  may  not  be  unique  in 
cases  where  m  =  n  was  taken  care  of  by  requiring  that  (a,b)  and  (b,a)  should  always 
be  taken  together  (in  the  critical  region). 

What  was  the  motivation  for  Barnard  to  include  this  third  condition?  Consider 
the  case  where  m  =  n  =  2,Xi  =  2,  and  x2  =  0.  It  follows  from  conditions  S  and 
C  alone  that  in  judging  the  significance  of  such  a  result  we  need  consider  only  the 
probability  of  this  result,  together  with  its  converse,  in  which  Xi  =  0  and  x2  =  2. 
In  this  case  the  attained  significance  level  is  given  by  292(1  -  9)2.  Suppose  we  had 
decided  beforehand  to  set  a  at  .05  .  Then,  if  .197  <  9  <  .803,  the  attained  significance 
level  would  be  greater  than  or  equal  to  a  and  hence  we  cannot  reject  Ho.  However, 
if  the  true  unknown  9  is  >  .803  or  <  .197,  the  attained  significance  level  would  be 
less  than  a  and  we  would  conclude  that  9\  ^  92.  There  is  a  unique  decision  only  if 
maxo<0<i  n{9)  is  less  than  a.  Therefore,  we  consider  only  those  critical  regions  that 
satisfy  maxo<0<i  n(9)  <  a. 

In  the  appendix  of  Barnard  (1947),  tables  of  attained  significance  levels  are  given 
for  the  sample  sizes  (7,7),  (8,6),  and  (9,5)  for  the  maximization  method.  Also,  for 
the  table  (7,7),  the  attained  significance  levels  corresponding  to  Fisher's  exact  test 
are  given.  For  each  observed  sample  point,  the  attained  significance  level  was  larger 
for  Fisher's  exact  test,  indicating  the  conservative  nature  of  this  test.  It  must  be 
noted  at  this  point  that  no  method  was  mentioned  of  how  to  solve  max0<e<i  tt(9). 
Indeed,  such  a  method  must  be  developed  since,  in  most  cases,  the  solution  is  not 
easily  obtained  (in  fact  cannot  be  obtained  except  for  very  small  samples)  using  the 
simple  method  of  derivatives.  There  was  not  much  research  done  in  this  area  in  the 
years  following  this  paper,  and  in  fact  Barnard  himself  rejected  the  method  in  favor 
of  Fisher's  exact  test  in  Barnard  (1949). 
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Boschloo  (1970)  suggested  a  method,  which  is  loosely  based  on  the  above  methods 
of  Barnard,  for  improving  the  power  of  Fisher's  exact  test  in  the  context  of  the  com- 
parative binomial  trial.  Again,  consider  Xx  and  X2  to  be  two  independent  binomial 
random  variables  with  parameters  (m,  9\)  and  (n,  92)  respectively.  Boschloo  considers 
both  the  one-sided  and  two-sided  cases  in  what  follows.  Let  the  attained  significance 
level  using  Fisher's  exact  test  be  denoted  by  cep,  and  the  corresponding  unconditional 
level  by  a{9),  i.e. 


where  t  is  the  total  number  of  successes,  H(xi,t)  is  as  defined  above,  C  is  the  set  of 
x\  values  over  which  the  attained  significance  level  is  to  be  calculated,  and  9  is  the 
common  value  of  9\  and  92  under  H0.  It  is  stated  that  numerical  calculations  show 


small  sample  sizes.  The  small  value  of  max^^i)  a{9),  the  size  of  the  test,  leads  to 
a  loss  of  power.  The  improvement  suggested  by  Boschloo  is  based  on  the  following 
idea.  Since  the  conditional  level  ap  results  in  a  small  unconditional  size,  a  raised 
conditional  level  will  be  used,  and  selected  in  such  a  way  that  max0e(O,i)  a(9)  <  a, 
the  desired  size  of  the  test.  The  raised  conditional  level  is  denoted  by  7.  7  is  made  as 
large  as  possible,  within  the  constraint  that  max^e(0,i)  a(9)  <  a.  The  calculations  of 
Fisher's  test  remain  unchanged,  only  the  level  of  significance  aF  is  raised  to  7,  where 


Using  a  raised  conditional  level  of  significance  in  most  cases  results  in  points  being 
added  to  the  critical  region,  which  in  turn  leads  to  an  increase  in  the  power  of  the 
test.  A  table  of  attained  power  is  given  for  the  common  sample  size  m  =  n  =  15, 
and  several  combinations  of  0X  and  92,  for  three  versions  of  Fisher's  exact  test.  These 
three  versions  are  1)  randomized,  2)  non-randomized,  and  3)  using  the  raised  level 


m+n-t 


that  max0e(O,i)  a(9)  is  often  smaller  that  | 


otp,  and  that  this  is  particularly  true  with 


7  =  sup{o!F  :  max  a(9)  <  a}. 
fle(o.i) 


8 


of  significance.  For  each  pair  of  6\  and  92  considered,  the  power  is  higher  for  3) 
than  for  2).  There  are  tables  of  7  values  for  sample  sizes  up  to  50,  and  the  a  values 
0.005,0.01,0.025,0.05,  and  0.10.  However  there  is  no  mention  of  the  technique  used 
to  determine  raaxee(0ii)         which  is  needed  to  find  7. 

McDonald  et  al.(1977)  considered  a  modification  of  Barnard's  maximization  method, 
using  as  the  critical  region  those  points  on  the  interior  of  the  corresponding  UMPU 
test's  critical  region  (i.e.  those  points  that  are  in  the  critical  region  but  are  not  used 
in  connection  with  randomization).  Again,  consider  two  independent  binomial  ran- 
dom variables  Xi  and  X2  with  parameters  (m,  6\)  and  (n,  62)  respectively.  For  testing 
H0  :  6\  =  02  versus  Ha  :  6\  <  92  ,  they  used  the  following  3  step  algorithm. 
Step  1  Determine  X*  such  that  F(X*)  =  H(x1  \  t)  <  a+,  but  F(X*+1)  >  a+, 
where  k  =  max(0,  t  -  n)  and 


H(x1  I  t)  = 


max(0,  t  —  n)  <  X\  <  min(£, m). 


This  step  is  used  to  form  a  critical  region,  based  on  a+,  of  the  form 

R  =  {(xi,t  -  xi)  :  xx  =  k,  ...,X*  and  t  =  0, 1, n  +  m}. 

Note  :  Although  it  is  not  formally  mentioned  in  the  paper,  we  believe  that  by  X* 
they  really  mean  X*(t).  It  does  not  make  sense  otherwise. 

Step  2  Numerically  determine  max7r(0)  =  a*,  where  tt(6)  is  the  null  power  function. 
They  do  not  describe  the  method  they  used  to  accomplish  this.  As  this  is  certainly 
not  a  trivial  computation,  the  method  should  have  been  given. 
Step  3  Iterate  on  values  of  a+  such  that  a*  <  a,  and  R  is  as  large  as  possible. 

The  motivation  for  this  paper  was  that  the  non-randomized  conditional  test  was 
too  conservative,  in  the  sense  that  the  actual  significance  level  can  be  one-forth  to 
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one-half  the  anticipated  value.  However,  no  comparison  of  the  results  of  the  above 
method  versus  the  exact  conditional  method  were  given. 

Haber  (1986)  proposed  a  one-sided  exact  test  based  on  the  unconditional  Z  statis- 
tic with  the  pooled  estimate  of  variance.  Let  0  <  Z\  <  z?  <  •  ■  •  <  Zk*  be  the 
different  positive  values  of  the  Z  statistic.  Define  Qk  =  maxg6(o,i)  Pr(Z  >  Zk).  For 
a  given  significance  level  a,  let  k*  be  an  integer  such  that  Qk*-i  >  a  and  Qk*  <  ct. 
Then,  the  exact  Z  test  is  defined  as  reject  H0  when  Z  >  c,  where  Zk*-\  <  c  <  z^- 
In  this  paper,  maxfle(0ii)  Pr(Z  >  Zk)  was  calculated  over  the  finite  set  of  points 
9  =  0.0, 0.01, 1.0  .  Power  comparisons  were  made  between  Fisher's  exact  test 
(non-randomized),  Fisher's  exact  test  (randomized),  and  the  exact  unconditional  Z- 
test  for  various  combinations  of  9\  and  92  with  two  combinations  of  m  and  n,  namely 
m  —  n  =  10  and  m  =  5,  n  =  15.  For  these  tables,  the  exact  Z-test  is  "considerably" 
more  powerful  than  Fisher's  exact  test  (non-randomized)  and  in  most  cases  its  power 
is  close  to  or  even  exceeds  Fisher's  exact  test  (randomized). 

Suissa  and  Shuster  (1985)  developed  a  rigorous  method  to  determine  the  maximum 
of  the  null  power  function  for  the  comparative  binomial  trial.  Their  methodology  was 
developed  to  calculate  sup/(0)  over  9  G  (0, 1),  where 

/(0)  =  £a^(l-0)Cl  (1.1) 

iec 

di,bi,  and  q  are  all  greater  than  or  equal  to  zero,  i  is  an  indexing  subscript  over  the 
whole  sample  space  S,  C  is  a  set  of  subscripts  for  which  the  related  sample  points 
belong  to  the  critical  region  defined  by  the  testing  procedure,  and  9  is  a  nuisance 
parameter. 

The  method  of  computing  sup/(0)  is  based  on  the  mean  value  theorem  of  differ- 
ential calculus  applied  to  successive  sub-intervals  of  the  unit  interval  (0,1).  As  will 
be  seen,  if  there  exists  a  bound  for  f'{9)  on  each  sub-interval,  we  can  use  the  mean 
value  theorem  to  get  an  upper  bound  for  f(9).  The  derivative  of  f{9)  is  given  by 
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f{6)  =  £  aAo^ii  -  er  -  aiCl9b>(i  -  ey>-1.  (1.2) 

iec 

This  is  a  linear  combination  of  binomial  terms  of  the  form  h{9)  =  9r(l  —  9)s~r. 
Thus,  for  any  given  interval  /  =  (a,  6),  where  0  <  a  <  6  <  1,  it  follows  that 


sup  h{9)  =  < 


h(a)  if  r-  <  a 
h(b)    if  ^  >  6 


(1.3) 


h{T-)  ifr-ei 

<        v  S  I  8 


and  also, 


vafh{9)  =  mm{h(a),h(b))  (1.4) 

An  upper  bound  for  f'(9)  can  then  be  obtained  by  substituting  the  right  hand 
side  of  (1.3)  for  each  positive  term  of  (1.2),  and  the  right  hand  side  of  (1.4)  for  each 
negative  term  of  (1.2).  A  lower  bound  for  f'(9)  can  be  obtained  on  (a,  b)  by  reversing 
these  substitutions.  Thus,  a  bound  M  for  |  f'(9)  \  on  (a,  b)  is  taken  as  the  larger  of 
the  two  bounds  in  absolute  value. 

We  can  now  obtain  an  upper  bound  for  f(9)  on  each  of  the  successive  intervals 
h  =  (0,0.01),/2  =  (0.01, 0.02),..., 7100  =  (0.99,1)  using  the  mean  value  theorem.  An 
upper  bound  for  f(9)  on  the  interval  (0, 1)  is  obtained  by  simply  using  the  maximum 
of  the  upper  bounds  from  the  sub-intervals.  For  each  =  1, 100  we  can  find  a 
Mj  such  that  |  f'(9)  \<  Mj  V9  G  /,.  By  the  mean  value  theorem  it  can  be  concluded 
that 

f(9)  G  (f(Sj)  -  0.005M,,  f(Sj)  +  0.005M;)  V0  e  Ij 

where  Sj  =  Thus,  for  each  Ij,  an  upper  bound  for  f(9),  9  €  Ij,  is  given  by 

fM  =  f(gj)  +0.005M,-,  so  that  f{9)  <  max^^ioo  fj(9)  .  Since  Mj  could  possibly 
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be  very  large  for  some  values  of  j,  the  above  bound  may  be  very  conservative.  This 
bound  can  be  improved  upon  to  produce  a  upper  bound  within  5  of  the  true  maximum 
in  the  following  manner.  Any  interval  Ij  for  which  fj(6)  >  max,/^)  +  8  must  be 
subdivided  into  2m>  sub-intervals  of  the  form  {Ijkj  :  kj  =  1, 2, , 2mJ }  where  rrij  is 
the  smallest  integer  such  that 

,  0.005Mjfc. 

f(s3kj)  +  — om         <    max(max/(si),  max  f{sjkj))  +  5 

Vfcj,  and  where  Sj^  —  ^  +  ^^^-i  is  the  midpoint  of  the  Ijkj,  and  Mjkj  is  the 
local  bound  for  |  f'(8)  |  in 

Using  the  methodology  developed  here,  it  is  possible  to  compute  the  size  of  any 
test  for  which  the  null  power  function  is  of  the  form  (1).  Let  X\  and  X2  be  indepen- 
dent binomial  random  variables  with  parameters  (n,  9\)  and  (n,  62)  respectively.  The 
problem  is  to  test,  at  level  a,  the  hypothesis  H0  :  61  =  62  versus  Ha  :  9\  <  92.  The 
methodology  they  developed  was  used  to  compute  the  size  of  the  test  based  on  the  Z 
statistic  with  unpooled  variance  estimator.  The  null  power  function  is  given  by 

=      £      (  r    )(  r    )  dXl+X2^  ~  Of"-*1-*2 
(xux2)ec  \Xl  J  \X2  ) 

where  9  is  the  nuisance  parameter  and  C  is  the  critical  region  defined  by  the  test 
statistic.  For  the  one  sided  test,  the  critical  region  defined  by  Z  is  given  by 

C  =  {{xi,  x2)  :  Z  >  z,  xi,  x2  =  0, 1, n  and  2  >  0}. 

For  an  oj  level  test,  the  critical  value  of  z,  call  it  2*,  satisfies  the  equation 

z*  =  inf{2  :  sup7r(#)  <  a). 

Now,  we  can  find  a*,  the  size  of  Z,  for  any  value  2.  Then,  for  a  test  of  significance 
level  a,  the  critical  value  2*  can  be  obtained  by  the  above  equation.  We  can  do  this  by 
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using  the  100(1  -  a)  percentile  point  of  the  standard  normal  as  a  starting  point.  This 
value  is  incremented  or  decremented  until  a*  <  a,  and  z*  is  taken  as  the  smallest 
value  which  satisfies  this  inequality.  This  procedure  was  done  using  a  FORTRAN 
program. 

For  n  =  10,...,  150  and  a  =  0.025  and  0.05,  z*  and  a*  (within  precision  0.001) 
were  computed  and  given.  Using  these  critical  values,  the  power  was  calculated  for 
Of  =  0.025  and  0.05  and  various  values  of  9X  and  92-  The  minimum  sample  size  required 
per  group  to  attain  a  power  of  1  —  (5  and  significance  level  a  can  be  computed  by 
solving 

n*  =  min{n  :  n(9u92)  >  1  -  /?} 

where  the  critical  region  is  based  on  2*,  a  function  of  n. 

Gail  and  Gart  (1973)  determined  the  sample  size  requirements  for  the  exact  con- 
ditional test.  There  are  also  two  other  commonly  used  methods.  The  first  is  based  on 
the  pooled  Z-test  and  the  second  is  based  on  the  variance  stabilizing  property  of  the 
arcsine  transformation  on  proportions.  It  was  shown  that  n*  agrees  quite  well  with 
the  sample  sizes  determined  by  the  pooled  Z  statistic  and  the  arcsine  method.  Also, 
n*  tends  to  be  smaller  than  the  sample  size  determined  using  the  exact  conditional 
test.  For  the  cases  considered  in  this  paper,  n  =  10(1)150  and  a  =  0.05,0.025,  and 
0.01,  the  exact  unconditional  Z-test  was  uniformly  more  powerful  than  Fisher's  exact 
test  (non-randomized).  A  similar  result  was  obtained  in  Suissa  and  Shuster  (1984), 
were  it  was  shown  that  the  exact  unconditional  Z-test  is  more  powerful  than  the  cor- 
responding UMPU  test  in  a  substantial  and  important  (meaning  that  the  values  of 
9i  and  92  were  not  near  zero  or  one)  region  of  the  alternative  parameter  space. 

A  rigorous  treatment  of  the  maximization  method  for  2  by  2  tables  in  the  context 
of  multinomial  sampling  is  given  in  Shuster  (1992).  For  the  matched  pairs  case  see 
Suissa  and  Shuster  (1991). 
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Berger  and  Boos  (1994)  proposed  a  method  which  could  possibly  simplify  calcu- 
lations and  increase  power.  In  the  presence  of  a  nuisance  parameter,  9,  they  redefine 
the  p-value  as 

Pfi  =  maxp(0)  +0 

where  Cp  is  a  1  —  /3  confidence  set  for  the  nuisance  parameter  when  H0  is  true.  If  the 
maximum  occurs  outside  the  confidence  set,  there  could  be  a  considerable  increase  in 
power.  It  is  claimed  that  this  adjusted  p-value  may  be  preferred  to  max0<=e  p{9),  where 
6  is  the  domain  of  the  nuisance  parameter  space,  since  it  can  be  computationally 
simpler,  and  it  is  also  sensible  to  determine  the  maximum  over  the  most  likely  regions 
of  0.  Some  examples  are  given,  most  notably  the  one  concerning  the  two-sided  test 
for  a  normal  mean  when  the  variance  is  unknown.  The  p-value,  as  a  function  of  a2,  is 
given  by  p{a2)  =  2<3>(-  |  2obs  |).  In  this  particular  case  we  cannot  use  max(T2>0  p(a2), 
since  it  will  always  equal  one.  Therefore,  in  this  particular  case,  some  modification 
of  the  maximization  method  was  necessary.  Also,  this  method  can  be  helpful  when 
comparing  two  independent  binomials,  when  one  sample  size  is  much  larger  than  the 
other. 

1.3    Motivation  for  an  Alternative  to  Conditional  Methods 

It  is  a  well  known  fact  that  the  non-randomized  version  of  Fisher's  exact  test 
tends  to  be  very  conservative  for  sparse  contingency  tables.  An  examination  of  the 
null  power  function  for  such  a  test  reveals  this.  Consider  a  2  by  2  contingency  table 
arising  from  the  comparative  binomial  trial,  where  the  sample  size  for  each  binomial 
population  is  10,  a  =  0.025,  and  the  hypothesis  to  be  tested  is  H0  :  9X  =  92  versus 
Ha  ■  9 1  ^  92.  The  graph  of  the  null  power  function  corresponding  to  this  testing 
scenario  is  given  in  Figure  1.1  (note:  all  Figures  appear  at  the  end  of  this  chapter). 
The  conservative  nature  of  Fisher's  exact  test  is  apparent.  In  fact,  the  maximum 
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power  attained  is  only  0.0127.  In  Figure  1.2,  the  graph  of  the  null  power  function  for 
the  common  sample  size  of  20  is  given.  Although  it  is  not  as  conservative  as  the  for 
the  previous  case,  the  attained  power  still  falls  well  short  of  a  =  0.025.  Indeed,  the 
maximum  power  is  only  0.0139. 

For  the  case  of  two  multinomial  populations,  we  are  interested  in  testing 
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at  a  —  0.025.  Here,  6\  =  (0n,  0i2, 0i3)'  and  02  =  (021,022  ,  023)'  are  the  cell  probabili- 
ties for  the  two  multinomial  populations.  Consider  a  testing  procedure  based  on  the 
non-randomized  version  of  Fisher's  exact  test,  where  the  number  of  trials  for  both 
populations  are  equal.  The  graphs  of  the  null  power  functions,  using  sample  sizes  of 
10  and  20,  are  given  in  Figures  1.3  and  1.4  respectively.  The  conservative  nature  of 
the  conditional  test  is  still  apparent.  However,  it  is  not  as  conservative  as  it  was  in 
the  2  by  2  context  (clearly  demonstrated  by  the  graphs). 

1.4    Statistical  Power  Comparisons 

Mehta  and  Hilton  (1993)  suggest  that  the  power  advantage  the  exact  uncondi- 
tional method  has  over  the  conditional  method  for  the  2  by  2  table  is  primarily  the 
result  of  the  discreteness  of  the  hypergeometric  distribution.  This  is  a  problem  which 
they  claim  would  not  be  a  factor  in  larger  contingency  tables.  In  fact,  they  hint  at 
the  possibility  that  a  conditional  test  could  be  more  powerful  than  a  corresponding 
unconditional  test.  However,  for  any  conditional  test,  there  is  always  an  uncondi- 
tional test  that  has  at  most  the  same  level  and  greater  than  or  equal  power.  This 
was  demonstrated  in  general  by  Suissa  and  Shuster  (1990)  and  Berger  (1994).  The 
argument  is  as  follows.  Let  pc  denote  the  conditional  p-value  of  the  test  statistic 
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employed  by  conditional  methods.  For  an  a  level  test,  the  null  hypothesis  is  rejected 
if  pc  <  a.  Now,  use  pc  as  the  test  statistic  for  the  exact  unconditional  procedure. 
Clearly,  Ph0(Pc  <!<*)<  a,  so  we  have  an  a  level  test.  However,  it  is  possible  that 
there  exists  a  a*  >  a,  such  that 

1.  Ph0(Pc  and 

2.  There  is  a  sample  point  with  a  <  pc  <  a*. 

In  this  case,  the  exact  unconditional  method  will  have  greater  power  at  level  a. 
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Figure  1.2  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  binomial  distributions  when  n  =  20. 
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Figure  1.3  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  10. 


Figure  1.4  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  20. 


CHAPTER  2 
MAXIMIZING  THE  NULL  POWER  FUNCTION 


2.1  Introduction 

In  this  chapter,  we  develop  a  method  for  determining  the  maximum,  with  respect 
to  0  —       92, Ok),  of  the  following  two  functions 

!■  /w  =  Enc(x)^,(i-^r"1"'  o<^<  1,^=0,1,-," 

A  t=l 

k  k  k 

2-  #)  =  EJIW.    O<0i<l,xi  =  O,l,..,n  ,  J>«  =  1,  J>  =  n 

where  C(x)  =  C(xi,a:2,  ...,a:jt)  and  D(x)  —  D(xi,x2:  ■■■,Xk)  are  positive  functions  of 
x  =  (xi,  x2, x/t)  ,  and  the  set  A  is  a  subset  of  the  set  {0, 1, n)  x  {0, 1, n)  x  •  ■ 
•  x  {0,1,  ...,n}  . 

The  function  f(0)  is  similar  to  the  null  power  function  for  the  comparison  of 
several  2  by  2  tables.  It  can  also  represent  the  likelihood  function  for  multivariate 
observations,  where  the  components  of  the  vector  are  independent  binomial  random 
variables,  and  hence  a  method  of  determining  this  functions  maximum  can  be  used  to 
find  maximum  likelihood  estimates  of  6  =  (6>i,  92, 0k).  The  function  g(6)  represents 
the  null  power  function  for  an  r  x  k  contingency  table  consisting  of  independent 
multinomial  random  variables.  The  size  of  the  test,  for  each  of  the  above  two  scenarios, 
is  given  by  supe6Q/  f(0)  and  supee09  g(0)  respectively,  where 

Sf  =  {(9ue2,...,9k)-0<91<1} 
20 
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and 

k 

eg  =  {(dl,92,...,eky,o<dl<i  $>  =  i} 

Note  that  since  f(0)  is  continuous  on  the  closed  set  0/  and  g(0)  is  continuous  on 
the  closed  set  Q0 


9' 


and 


sup/(0)=max/(0) 


supg(0)  =  maxg(0). 

0eQ„  "e09 


The  remaining  parts  of  this  chapter  deal  with  the  development  of  methodology  to 
determine  maxg€e/  f(0)  and  max^6e9  g(0)-  This  technique  will  be  used  in  two  cases. 
In  chapter  3  we  consider  the  2  by  3  table,  and  in  chapter  4  the  case  of  two  2  by  2 
tables,  i.e.  the  2x2x2  case. 


2.2    Finding  the  Maximum  of  f(0) 

Before  beginning  we  introduce  the  following  notation  to  make  the  presentation 
easier  to  read.  Let 

a.  Bk  =  <g>£=1  B,  where  B  is  any  set, 

b.  =  {0, 1, 2,  ...,n-l}*, 

c-  Ijij2,-,jk,ns  =  ®?=i  ^  where  j%  G  Ni,n«-i  and  s  and  n  are  positive  integers. 
Also,  for  each  set  /j1ja,...j4,n«  we  have  the  corresponding  sets 

d-  ?tlit2>...itjkjS)n,0  =  [^j^,  where  t<  =  n^s0)  -,n(ji  +  l)s0  -  1  and  s0  is  a 
positive  integer. 
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Lemma  2.1  V(ji,  j2, jfc)  e  Nk,na-i, 


k 

E 

y4  »=1 


A  i=l 


where  0  G  Ijltj2,...,jk,n8,  and 


Ts  ( Jt  >  t ) 


2L  >  £i 
ns  —  n 


7i  +  l  JL  <f  £ 

ns 


< 


pf.  Let  (ji,j2,---,jk)  G  Nk,ns-i  •  Since  C(x)0tXi(l  -  0j)n  x«  has  a  unique  maximum 
at  Qi  —  ^  (when  ^  G  [0, 1])  and  is  non-negative,  we  have  that 


Vi,    C(x)^(l-^)n_x<  <C7(x)7i(j<,^r(l-7f(j4la5i))*^    (*lf*a,  ...,*»)  €  i4 


i=l  i=l 

=>  Encw^(i-^^<En^w7.(»,^ni-7.(»i,«i))^ 

4  i=l  A  i=l 

Since  the  choice  of  (ju  j2, jr*)  was  arbitrary,  the  result  holds 

V(?i,ia,...,ijfc)  €A»,«,_i.  □ 

Using  Lemma  2.1,  we  have  an  upper  bound  for  on  any  rectangle  Ijuj2,...,jk,ns 
or  Th,j2,-,jk,s,ns0-  Let  UI(juj2,-;jk;s)  and  UT(ti,t2,  ...,tk;  s,  s0)  denote  the  upper 
bounds  of  /(0)  on  the  rectangles  Ijuj2,...,Jk,ns  and  rtl)t2,...jtjki,ifMo)  respectively.  Let 
M/,s  denote  the  maximum  of  the  UI(ji,j2,>~,jk;s),8  over  the  indices  (Juh,—,jk), 
and  let  MTtStSo  denote  the  maximum  of  the  tfr(*i>  *2> —>**;*»  *o)'s  over  the  indices 
(ti,*2,  Since  \J Ij1ji,...ji,tu  =  ©/>  we  have  that 

<  M/>a       V  0  €  6/  and  s  €  N+ 
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where  iV+  is  the  set  of  positive  integers.  Therefore  MI<S  is  an  upper  bound  for  f(6) 
on  [0,  l]k  .  However,  it  may  be  the  case  that  MIiS  is  not  very  tight,  i.e. 
|  MIiS  -  maxee&f  f(0)  |  is  relatively  large. 

Suppose  it  is  desired  to  have  an  upper  bound  for  f(0)  within  5  of  the  true  maxi- 
mum. We  can  do  this  as  follows.  Consider  any  rectangle  Ijuj2,...jk,ns  for  which 

Ui(ji,32,-,jk;s)>    max    f(xi)  +  8  (2.1) 

l<i<(ns)k 

where     is  the  midpoint  of  Ij1,j2,...,jk,ns  and  i  is  an  index  for  the  sets  {/jlvj2,— j*,n* 

}.  This 

rectangle  must  be  further  subdivided  into  the  sets  Ttlj2,...,tk,s,ns0-  F°r  each  rectangle, 
Ijl,j2,...,jk,nsi  which  satisfies  (2.1),  we  want  to  find  the  smallest  positive  integer  s0  such 
that 

Mt,s,s0  <  max  <^    max    /(x^),    max  /(xj)>+6 

U<i<(ns)*  l<l<(ns0)k  J 

where  x;  is  the  midpoint  of  Ttut2,...,th,s,na0,  and  /  is  an  index  for  the  sets  {7flit2i...,tfciS)nso}. 
Then,  we  can  determine  an  upper  bound  for  f(0)  that  is  within  S  of  the  true  maxi- 
mum. Namely, 

max{MT,S)S01,  MT^S02, MT<SiSOg,  M) 

where  g  is  the  number  of  rectangles  Ijuj2,...,jk,ns,  for  which  (2.1)  holds,  and  M  is  the 
maximum  of  the  Ui(ji,j2,---,jk',sys  over  the  set  of  indices  (Ji,j2,—,jk)  for  which 
(2.1)  does  not  hold. 

It  remains  to  show  that  it  is  indeed  possible  to  find  an  s0  for  each  rectangle, 
Ijij2,-jk,ns,  where  UI(j1,j2,—,jk,s)  >  max,  /(x,)  +  5.  We  will  now  demonstrate 
that  there  exists  a  subsequence  of  positive  integers,  {si,  s2, ...},  such  that  M/)Sl,  M/iS2, 
M/jS3,  ...  converges  to  maxe€e/  f(9),  and  hence  s0  exists. 
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Theorem  2.1  If  S  =  {si,s2,...}  is  a  sequence  of  positive  integers  such  that  Sj  is  a 
multiple  of    Vz  <  j,  then 

lim  Mi.  =  max  f(0) 
i-*oo    ',s'     flee, J  v  ' 

pf.  Suppose  0  €  Ijida,...jkW  an<^      -^/O'l)  J2)  — »i*>  s)  be  defined  as  follows 

fc 


A  i=l 


where 


'  Jt  +  l  ii_  >  £i 

ns      '      ns  —  n 


1L 
n.s 


ns  n 


Similar  to  Lemma  2.1,  we  have  that  Lj{ju  j2, jk;  s)  <  f(0),  V0  G  Ijuj2,...,jk,ns,  and 
hence 

L/)S<max/(0)<M/jS  (2.2) 

where  LIyS  is  the  maximum  of  the  LI(ji,j2,—,jk',s),s  over  the  indices  (ji,j2,  j*). 
Define  the  following 

«*  =  (7«(ji,^))X<  (1  -  Jstiux^y-^  and  6,  =  frfaxi))*  (1  - 
4  =  nUi  C(x)am  and      =  UL=i  C(*)bm 
A0  =  B0  =  1. 

Now,  select  any  point  such  that  s  6  5,  and  use  it  as  a  "reference 

point".  We  can  write  t//(j1}  j2,        a)  -  L/0"i,  j2) j*;s)  as 


.4  k  i=i 
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and,  note  that 


CLi-bi=  (%(ji,Xi))Xi(l  -  js(ji,Xi))n  Xt 


1  - 


'Ss(ji,Xi)\Xi  ( 1  -  Ss(ji,XiY 


It  must  be  the  case  that  either 


Js(ji,Xi)J  \l-%{ji,xi)i 


Ji    \    f    ns-ji  \ 


Ji  +  l)  \ns-ji-li 


or, 


r8,{ji,Xi)\Xi  ( 1  -  6a{jj,XiY 

JaUuXi))      VI  -%(Ji^Xi)i 


'ji  +  l\Xi  ( ns-jj  -  1 
,  ji   J    V   ns-ji  , 


Vi  =  1, 2, A;.  For  an  arbitrary  /,  at  our  reference  point,  but  with  the  partition  made 
finer  (by  using  st  in  place  of  s,  where  si  E  S  and  s/  >  s)  we  have  that 


Ji    \    I    ns-ji  \ 


Ji  +  l  J  \ns-ji-l 


becomes 


\  %i  /         9       •         \  n—Xi 

jiSi    \    /    nsz  -  jiSi  \ 


JiSi  +  sl  \ns2-jlsi-s 


and, 


.    ji    J    I  ns-ji 


becomes 


f  .  \  Xi  /      o        •  \  n—Xi 

jiSi  +  s\      nsz  -  jiSt  -  s\ 


JiSi 


ns*  -  jiSt 


Clearly, 

_jisi_\  (  nf-m  \  =1 

3iSi  +  sJ  \ns2-jiSi-sJ 

and, 

Jjsi  +  A      ns  -3i8l-s\  =1 

jiSi    J  \    ns2  -  jiS[  I 
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since  Si  -4  oo  as  /  — >  oo.  So,  for  any  fixed  point  (£,  we  choose, 

£//(T ,  *f ,     *f I  *i)  -  >  T ,  *i)  -+  0  as  I  ->  oo  .  It  follows  that 

Ve  >  0,  3m  €  TV  such  that 

Ui(ji,j2,-,jk;si)  -  Li(jl,j2,-,jk;st)  <  e    V(ji,j2,-,jk)  e  iV^^-i 

whenever  /  >  m.  This  implies  that  M/  S;  —  — >■  0,  as  /  — >•  oo,  which  proves  the 
claim.  □ 


2.3    Finding  the  Maximum  of  g(0) 

The  domain  of  g(0)  is  different  from  that  of  f(0).  As  a  result,  we  need  to  "adjust" 
our  notation  as  follows 

a-  r n,n,- jk,ns  =  {®L  [i.^lJn©*  where  jx  €  NliM-X  and  ji  =  0,...,ns  -  1  - 
S'=l  ii>  and  for  each  set  I'juj2l-,jk,ns  we  have  the  corresponding  sets 
b-  T'hJ2,..Jlc,s,ns0  =  {®f=i  [jfe,  ^S] }  D  09  where    =  njiSo, n(j,  +  l)s0  -  1  Vi 
c-  U'i(juh,-,jk,s),  U'T(j1,j2,-,jk,s),  M'Ita,M'T,8,8o  are  analogous  to 
Ui(h,32,-,jk,s),  UT(jiJ2,-,jk,s),  and  M/iS,MTiSiSo,  but  are  with  respect  to  the 

qpfo    V  .  onrl  T" 

For  what  follows,  it  is  important  to  realize  that  the  function  t(0)  =  C(x)  []*=„  0iXi, 
with  the  constraints  £(9,  =  1  and  E^  =  n  is  maximized  at  §i  =  Since  this  fact 
is  not  so  obvious,  it  is  demonstrated  using  Lagrange  multipliers  as  follows. 

Let  h{0)  =  \ogt{0)  =  logC(x)  +  Eti  Xi  \og6i  and  1(0)  =  E  ^  -  1  =  0.  Now, 
V/i(0)  =  AVZ(0) 
=►  £  ^ei  =  A  £  e* 
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-±  £l  —  22.  —  .  .  .  —  *k 
Bi  ~  62  ~  9k 

=>0i  =  0k%   V«  =  l,...,*-1 

^  ek  =  since  £0;  =  1  and  £x8  =  n.  Thus,  i(0)  is  maximized  at  Ot  =  *  for 
i  =  1,2,..., k.  This  is  indeed  a  global  maximum,  since  either  2(0,0,  ...,0)  =  0  or 
i(l,  1, 1)  =  0,  and  t(0)  is  non-negative. 

Lemma  2.2  Let  0  G  I'juj2,...,jk,ns-  Then 

Jt  A: 

/l     1=1  y4     1  =  1 

where 


7s  (ji) 


ns 

Ji  +  1 
ns 


JL  > 
ns  — 


ii.  <  Si 


pf.  similar  to  Lemma  2.1. 

Similar  to  section  2.2,  we  can  use  Lemma  2.2  and  the  fact  that  U  I' juj2,...,jk,ns  —  Qg, 
to  get  an  upper  bound  for  g(6)  on  Qg.  Namely 

g(0)  <  m'/,s    v  e  e  eg ,  s  g  n+. 

Again,  it  may  be  the  case  that  this  upper  bound  is  not  very  tight.  We  can  get  an  upper 
bound  within  S  of  the  true  maximum  so  long  as  there  exists  a  sequence  of  positive 
integers,  {si,s2,...},  such  that  M'IiSl,  M'I<S2,  M'I<S3,  ...  converges  to  maxflee9  g{0). 

Theorem  2.2  If  5  =  {si,s2,...}  is  a  sequence  of  positive  integers  such  that  Sj  is  a 
multiple  of  Si  V«  <  j,  then 

lim  M'r  Si  =  maxq(0) 
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pf.  Suppose  0  G  I'jltj,,...jk,ns,  and  let  L't(jl,j2,  -Jk\s)  be  denned  as  follows 


A  i=l 


where 


'  J'i+l  ji.  >  Si 

ns      '      ns  —  n 


&s  {ji )  ) 


A. 

IIS 


ns  n 


Similar  to  lemma  2.2,  we  have  that  L'itt(ji,j2,...,jk;s)  <  ^  e  ^'iij2.  -j*.n*> 

and  hence 

///,,  <  max        <  M'/)S 

where  Z/s  is  the  maximum  of  the  L'(ji,j2, jfc;s)'s  over  the  indices  (ji,j2,  ■  ■■,3k)- 
Define  the  following 

a.  =  (T.O'i^i))^  and  bi  =  ($s(ji,Xi)T' 

A  =  nUi  C(x)am  and  ^  =  nUi  C{x)bm 

A0  =  BQ  =  1 

Select  any  point  ii, ... ?  it)  such  that  s  e  5,  and  use  it  as  a  "reference  point".  We 
can  write  U'I(Jl,j2,...,jk;s)  -  L'I(j1,j2,...Jk;s)  as 


and,  note  that 


«t  -  6i  =  (7*0«^t))X< 


'$s{ji,Xi)\Xi 

Js(ji,Xi)l 


29 


It  must  be  the  case  that  either 


$s  {ji  j  %i )  \  (  ji 


\%(ji,Xi)J  \ji  +  1J 

or, 

fs.Ui,xi)y  =  (ji±i)Xi 

\%(ji,Xi))  \     ji  J 

Vi  =  1, 2, k.  For  an  arbitrary  /,  at  our  reference  point,  but  with  the  partition  made 
finer  (by  using  S[  in  place  of  s,  where  s;  €  S  and  sj  >  s)  we  have  that 


becomes 


and, 


^-r—  ]       becomes     (  ^',St    S  ) 

<    Ji    J  V    Jt*l  / 


Clearly, 


and, 


lim  f-^r=i 


lim  i'^r=i 


since  s,  -»  oo  as  /  ->  oo.  So,  for  any  fixed  point  l^,      ....      we  choose, 
U'i(i&,  *f, Sl)  _  ^(Aft,  fet, Aa;  5|)       o  as  /  ->  oo  .  Thus,  M',,,,  - 
— >  0  as  /  — >  oo,  which  proves  the  claim.  □ 


CHAPTER  3 
THE  2x3  CONTINGENCY  TABLE 


3.1  Introduction 

In  chapter  1  a  methodology  for  determining  the  maximum  of  the  function 

5>0a'(l-0)6' 

R 

with  respect  to  9,  as  developed  by  Suissa  and  Shuster  (1985),  was  given.  In  the 
papers  Suissa  and  Shuster  (1985)  and  Suissa  and  Shuster  (1991),  this  method  was 
used  to  determine  exact  unconditional  critical  regions  for  hypothesis  tests  involving 
independent  and  correlated  proportions.  These  critical  regions  can  then  be  used  to 
determine  the  sample  size  requirements  necessary  to  achieve  a  given  power. 

The  results  were  favorable,  as  each  exact  unconditional  test  was  uniformly  more 
powerful  than  its  conditional  counterpart  for  all  of  the  cases  considered.  Also,  the 
sample  size  requirements  to  achieve  a  given  power  were  almost  always  smaller  using 
exact  unconditional  methods  as  compared  with  conditional  ones.  The  aim  of  the  cur- 
rent chapter  is  to  develop  exact  unconditional  critical  regions,  using  the  methodology 
developed  in  chapter  2,  and  determine  sample  size  requirements  for  the  comparison  of 
two  independent  trinomial  populations.  We  begin  with  a  description  of  the  problem. 

Let  X  and  Y  be  independent  trinomial  random  variables  with  parameters  (n,  0X) 
and  (n,0y)  respectively,  where  0X  =  (#n,  #i2,  M  and  0y  =  (6>2i,  022,  M  •  The  equal 
sample  size  case  is  considered  here,  because  Lachin  (1977)  has  shown  this  is  optimal, 
for  a  given  total  sample  size.  We  are  interested  in  determining  whether  0X  =  0y, 
i.e.  the  homogeneity  of  the  two  trinomial  populations.  The  hypothesis  of  interest  is 
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formally  stated  as 

H0  :0x  =  0y  vs.  Ha:0x^6y 

and  is  to  be  tested  at  a  given  level  a.  The  outcome  of  such  an  experiment  can  be 
represented  in  the  following  2x3  contingency  table 


■'■■2 

yi 

2/2 

2/3 

The  probability  of  observing  this  table  is  given  by 


where  Q{  is  the  common  value  of  $u  and  02t  under  H0,  for  i  =  1, 2, 3.  Since  £  @i  =  1, 
£  Xj  =  n,  and  £  Ui  =  ft>  we  can  write  this  probability  as 

P(x,y)  ^  (  V      n      1  IT  ^"(l  -  0X  -  02)2"-*'-^-^ 

a  function  of  the  two  nuisance  parameters  #x  and  #2-  Because  of  the  nuisance  pa- 
rameters, asymptotic  or  conditional  methods  have  traditionally  been  used  to  perform 
the  above  hypothesis  test.  There  have  been  no  attempts  made  to  compute  the  exact 
unconditional  test  size  because  of  the  computational  complexity  involved.  Hence,  the 
maximization  method  of  Barnard  (1947)  has  not  been  used  for  this  problem.  In  this 
chapter,  methodology  developed  in  chapter  2  will  be  used  to  compute  exact  uncondi- 
tional test  sizes,  allowing  us  to  form  critical  regions  using  the  maximization  method. 
First,  we  briefly  review  the  conditional  and  asymptotic  approaches  to  this  problem. 
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3.2    The  Conditional  Approach  to  the  2  x  3  Contingency  Table 

We  start  with  the  conditional  method,  which  is  usually  referred  to  as  Fisher's 
Exact  Test  (FET).  In  this  method,  conditioning  is  used  to  completely  specify  the 
relevant  probability  distribution  when  H0  is  true.  Since  the  trinomial  probability 
distribution  is  a  member  of  the  exponential  family,  the  nuisance  parameters  can  be 
eliminated  by  conditioning  on  appropriate  sufficient  statistics.  The  joint  distribution 
of  the  two  trinomial  random  variables  X  and  Y  is  given  by 


P(x,y)  = 


This  probability  distribution  will  be  completely  specified  under  H0,  i.e.  no  nui- 
sance parameters,  if  we  condition  on  X\  +  Yx  and  X2  +  Y2  .  Before  deriving  the 
conditional  distribution,  define  the  following 
t  =  (ti,t2)  ,  where  £,  =  Xi  +  y{  for  i  =  1, 2 

At  =  {(x,y) :  Xi  +  yi  =  <<  i  =  l,2  Xi,y,->0  £x,  =  n  Eyi  =  n} 

ai  =  max(0,  ti  —  n) 

a2  =  min(n,  t\) 

bi  =  max(0,  tx  +  t2  -  n  -  xx) 

b2  —  min(n  —  Xi,t2). 
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Then, 

P(X  =  x|Xi  +  Yi  =  tuX2  +  Y2  =  t2) 


D(x,  y)  exp  {eLi  U  log  0*  +  ZU  x>  log  (fe)  } 
J2At  D(x>,  y')  exp  {E?=1    log  0*  +  £?=i  xj  log  (  £})  } 

D(x,y)exp{ELi^log(fe)} 
E4t^(x',y')exp{EL^log(fc)} 


n  \  (  n  —  Xi  \  (  n  \  I  n  —  t\  +  X\ 
x\  I  \     x2     /  V  *i  —  2/i  /  V     *2  -  ^2 


exp{EL^log(^)} 


/  nWn-xiW      n      \  /  n-ti+ii  j 
V      /  \     ^2     /  V  h  -  xx  J  \     t2-x2  J 
(  n  \  (  n  -  xi  \  (      n      \  (  n  -  tx  +  xr' 
^At  Ui'     \  )  Wi  -  an'  j  I     *2  -  x2' 


I  n  \  /  n  -  xi  \  /      n      \  (  n  -  tx  +  xx  \ 
va2      V62      (  n  \  f  n-X!1  \  (      n      \  f  n  -  U  +  Xi 


Finally,  since 


n  \  (  n  -  xi  \  (      n      \  (  n  -  ti  +  X\  \      (  2n\  (  2n-t\ 


we  have  that 


n  \  (  n  —  Xi  \  f      n      \  (  n  —  ti  +  X\ 


P(X  =  x\Xi+Yi  =  tu  X2+Y2  =  t2)  ^ 


X2       I  \  ti  -  Xi   I  \      t2  -  x2 


2n\(2n-tx 
h  I  \  t2 
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where 

max(0,  U  -n)  <X\<  min(n,  t\)  and, 
max(0,fi  +t2  -  n  -  Xi)  <  x2  <  min(n  -  Xi,t2). 

We  can  now  determine  the  critical  region  for  a  test  of  level  a.   Let  x*  be  an 
arbitrary  value  of  x  =  (xi,x2,x3)  and  denote  the  conditional  probability 
P(X  =  x|Xi  +  Yi  =  U,X2  +  y2  =  *2)  under  tf0  by  P0(x|t).  Define 

Po(x*|t)  =  £P0(x|t) 

Bt 

where 

Pt  =  [(x,y):P0(x|t)<P0(x*|t)]  . 
For  each  value  of  t  and  n,  let 

C„,t(a)  =  [(x)y):p0(x|t)<a]  . 

Then,  from  Freeman  and  Halton  (1951),  the  critical  region  is  given  by 

Rn(ot)  =  \JCn,t{a)  . 
t 

If  the  size  of  the  test  is  strictly  less  than  a,  which  will  usually  be  the  case  due 
to  the  discrete  null  distribution,  a  randomization  can  be  used  to  achieve  size  a.  The 
randomized  test  is  given  by 


0(x|t)  =  < 


fl  xeCn>t(a) 
7t   x  =  xt 
0  otherwise 


where  xt  e  At\Cnit{a)  and  po(xt|t)  <  po(x|t)  Vx  e  At\Cntt(a).  Also, 

a  ~  EcB,tW  ^o(x|t) 


7t 


Po(xt|t) 


However,  randomized  tests  are  not  typically  used.  Without  the  randomization,  the 
conditional  test  is  conservative,  which  was  seen  in  the  examples  from  chapter  1. 
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3.3    Determining  Sample  Sizes  Using  the  Conditional  Approach 

Two  steps  are  necessary  in  order  to  determine  the  sample  size  required  to  achieve  a 
given  power.  First,  for  every  value  of  n,  the  common  sample  size  of  the  two  trinomial 
populations,  critical  regions  are  determined  corresponding  to  the  level  of  the  test. 
This  is  Rn(a),  which  was  given  in  the  previous  section.  Second,  for  parameter  values 
under  the  alternative,  the  power  is  calculated  iteratively  with  respect  to  n  until  the 
required  power  is  attained. 

However,  the  critical  regions  get  quite  large  with  increasing  n.  This  makes  their 
determination  formidable,  since  the  large  number  of  calculations  requires  too  much 
computing  time.  To  see  this,  note  that  for  a  given  value  of  n,  the  possible  number  of 
2x3  contingency  tables  is 

'(n  +  l)(n  +  2)l2 
2 

which  increases  rapidly  with  n. 

A  simplification  can  be  made  which  reduces  the  calculations  and  amount  of  com- 
puter time  considerably.  In  what  follows,  it  will  be  shown  that  it  is  possible  to  simplify 
the  calculation  of  the  critical  region  by  condensing  the  sample  space.  This  will  sig- 
nificantly reduce  computing  time,  which  is  very  helpful  when  n  gets  moderately  large 
(>  25).  Define  the  following  set: 

Un  —  the  set  of  all  2  x  3  tables  such  that  any  table  in  Un  is  not  a  row  or  column 

permutation  of  any  other  table  in  Un,  and  each  table  is  ordered  by  column  totals 

and  first  row  elements. 
Although  the  description  of  the  set  Un  is  relatively  simple  to  understand,  it  is  crucial 
that  it  is  understood  correctly  because  of  it's  relevance  in  what  follows.  Therefore,  a 
example  is  given  in  Appendix  A  for  the  case  n  =  3,  i.e.  U3.  A  Fortran  program  was 
written  which  determines  Un  for  any  value  of  n,  and  also  determines  the  number  of 
unique  tables,  henceforth  referred  to  as  representations,  each  element  of  Un  generates. 
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For  a  given  value  of  n,  Un  is  between  ^  and  ^  the  size  of  the  original  sample 
space.  Sample  size  determinations  can  be  made  using  only  the  tables  in  Un,  which  is 
now  demonstrated. 

When  calculating  the  conditional  p- value  for  a  given  table  under  H0,  considering 
the  entire  set  of  tables,  not  just  those  belonging  to  Un,  the  usual  method  would  be 
to  first  calculate  the  conditional  probability  of  the  given  table, 

n\fn  —  xi\f     n     \  f  n  -  U  +  x\  \ 
Xi  )  \     x2     )  \  h  -  Xi  )  \     t2-x2  J 
(  2n  \  (  2n-U\ 

U/l    **  ) 
which  we  denoted  as  P0(x|t).  Then,  the  conditional  p- value  is  given  by 

Po((x',y')G  At :  P0(x'|t)  <  P0(x|t)). 

If  this  p- value  is  less  then  a  specified  a,  the  point  (x,  y)  will  be  included  in  the  critical 
region. 

However,  when  using  only  those  tables  in  Un,  a  modification  of  Po(x|t)  is  required 
when  calculating  the  p-value.  This  change  is  necessary  since  a  given  table  in  Un 
will  typically  have  many  representations.  Therefore,  we  must  determine  appropriate 
probabilities  associated  with  tables  in  Un,  which  will  be  denoted  P0*(x|t).  To  calculate 
P0*(x|t),  we  partition  Un  into  the  following  three  groups  of  tables: 

3.1.  t\  =  t%  =  £3, 

3.2.  ti  =  tj  for  some  i  /  j  and  3.1  is  false,  or 

3.3.  U  ^  tj  for  i,  j  =  1, 2, 3  and  i  /  j 

and  consider  each  group  separately.  First,  consider  case  3.1,  where  3  "types"  of  tables 
generate  ^1=^2  =  ^3: 

3.1a.  All  elements  in  the  table  are  the  same,  for  example 
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In  this  case  F0*(x|t)  =  P0(x|t). 

3.1b.  3i  €  {1, 2, 3}  such  that  x,  =  yu  but  Xj  +  yj  Vj  G  {1,2, 3}/i.  For  example 


8 

5 

2 

2 

5 

8 

In  this  case  P0*(x|t)  =  6  x  P0(x|t). 

3.1c.  Not  cases  3.1a  or  3.1b  .  For  example, 


9 

4 

2 

1 

6 

8 

In  this  case  P0* (x|t)  =  12  x  P0(x|t). 

This  exhausts  the  possibilities  for  case  1.1.  For  case  1.2,  there  are  four  "types"  of 
tables: 

3.2a.  3i  €  {1,2,3}  such  that  x{  =  and  U  ^  tj  for  i  /  j  and  j  G  {l,2,3}/i  For 
example, 


CO 

2 

3 

2 

8 

3 

In  this  case  P0*(x|t)  =  2x  P0(x|t). 
3.2b.  X2  =  £3  and  y2  =  y3-  For  example, 


10 

0 

0 

0 

5 

5 

In  this  case  P0* (x|t)  =  2  x  P0(x|t). 
3.2c.  x\  =  x2  and  yi  =  y2.  For  example, 


10 

10 

0 

4 

4 

12 
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In  this  case  P„*(x|t)  =  2  x  P0(x|t). 

3.2d.  Not  3.2a,  3.2b,  or  3.2c  .  For  example, 


9 

2 

3 

1 

8 

5 

In  this  case  P0*(x|t)  =  4  x  P0(x|t). 

For  case  3.3  there  are  just  two  "types"  of  tables: 

3.3a.  Xi  ^  i/i,  for  i  =  1,2,3.  For  example, 


In  this  case  P„*(x|t)  =  2  x  P0(x|t). 


10 

0 

1 

3 

8 

0 

mple, 

8 

5 

1 

8 

5 

1 

In  this  case  P0*(x|t)  —  A(x|t). 

For  cases  3.2  and  3.3  it  is  crucial  that  the  tables  be  ordered  as  t\  >  t2  >  t3  and 
ti  >  t2  >  tz  respectively,  in  order  that  what  follows  is  correct. 

The  p-value  for  any  table  in  Un  with  entries  xx,  x2,  X3,  yi,  y<i,  and  2/3  is  calculated 


as 


where, 


P0((x/,y/)€  A*:n*(x'|t)<P0*Wt)) 


A;  =  Unf\At 


The  critical  region  will  contain  all  those  tables  in  Un  that  have  corresponding  p- values 
<  «. 

Now  that  the  critical  region  can  be  formed  for  a  given  value  of  a,  it  is  possible 
to  make  power  calculations.  Before  presenting  the  form  of  the  power  function,  define 
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the  following 

cx{ex,9v)  =  «G«3«M(i  - en  - e12)2n-x<-*i(i  - e2l  -  en?n-"-*. 
c2(9x,0y)  =  ffiflg^flgfi  -  oxx  -  el2)2n~y^(i  -  e21  -  e22)2n~x^. 

Cid  =  {(1,2),  (2,1),  (2, 3),  (3, 2),  (1,3),  (3,1)} 

A,,  =  {(1,2),  (2,1),  (3, 2)} 
£^  =  {(1,2),  (3, 2),  (1,3)} 

For  power  calculations,  if  the  table  has  6  representations,  then  its  contribution  will 
be  one  of  the  following 

1.  For  contingency  tables  of  the  form  3.1fc,  3.2a,  and  3.36,  the  contribution  is 

2.  For  contingency  tables  of  the  form  3.26  the  contribution  is 


D(x,y) 


/]  C\(0X,  0y)  +  ^2  C2{0X,  By) 


3.  For  contingency  tables  of  the  form  3.2c  the  contribution  is 


£>(x,y) 


^  C\{BX,  By)  +  22  C2(BX,  By) 


and,  if  the  table  has  12  representations,  then  its  contribution  to  the  power  will  be 


£(x,y) 


52C1(Bx,By)  +  J2C2(Bx,By) 


Now,  proceed  by  iterations  on  n  as  follows.  For  given  values  of  Bx  and  By  in  the 
alternative  parameter  space  Ha,  calculate  the  power  iteratively,  starting  at  the  lowest 
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value  of  n  (=  10),  increasing  n  by  one  unit  until  we  achieve  the  desired  power.  This 
is  the  required  sample  size,  which  can  be  expressed  as 

nc  =  inf{n  :  7r„(0x,  0y)  >  ($} 

where  the  subscript  c  denotes  that  the  sample  size  is  based  on  conditional  methods, 
(3  is  the  desired  power,  and  7rn(0x,  0y)  represents  the  power  function  for  a  given  n. 

Before  beginning  the  section  on  asymptotic  methods,  two  important  questions 
must  be  answered  to  ensure  the  validity  of  the  above  procedure. 

1.  When  calculating  P0(x|t),  why  can  several  permutations  of  the  same  table  be 
accounted  for  multiplicatively? 

2.  When  determining  which  tables  from  Un  are  to  be  included  in  the  critical  region, 
tables  such  as 


0 

10 

0 

2 

1 

7 

are  not  considered.  However,  if  the  table 


10 

0 

0 

1 

7 

2 

is  in  the  critical  region,  then  the  former  table  will  be  used  in  the  power  calculation. 
Why  is  this? 

The  answer  to  both  of  these  questions  is  that  for  a  given  table,  conditional  proba- 
bilities are  invariant  to  row  and  column  permutations.  Consider  a  table  with  entries 

x2,  X3,  yi,  2/2,  and  y3.  It  is  clear  that  there  is  invariance  with  respect  to  a  row 
permutation,  since 

n  \  f  n  -  xi  \  (     n  \(n-ti+xi\ 
*i  )  \    x2    jyti-Xijy    h-x2  J 
2n  \  (  2n  -  tx  \ 
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obviously  equals 

{  n  \  f  n-yi\  /  n  \  (  n  -  h  +  yx  \ 
\  Vx  J  V    Vi     J  \ti-Vi  J  \    t2-y2  J 

V  *i  A  ) 

where  Xi  +  y\  =  *i  and  x2  +  y2  =  t2.  However,  the  fact  that  this  is  true  for  column 
permutations  is  not  so  obvious.  Consider  the  tables 


and 


Xi 

X2 

£3 

yi 

Vl 

yz 

Xi 

m 

2/3 

V2 

with  conditional  probabilities  given  by 


n 

Xi 


n  —  X\ 
x2 


n 

t\  -  Xi 


n  —  t\  +  X\ 

h  ~  x2 


2nU  2n  -  ti 

h  )  \  t2 


and 


fn\ln-xi\l      n      \  I  n  -ti  +  xi  \ 

V  Xi  A      X3      A       ~  Xl  A      *3  ~  ^3  J 

/  2nW  2n  -  \ 

1*1/1       *3  J 

respectively.  After  some  algebra,  the  ratio  of  the  above  quantities  reduces  to 

x3\  (n-xi-  x3)\  (t3  -  x3)\  (n  -  tj,  -  tz  +  xx  +  x3)\ t2\  (2n  -     -  f2)! 
x2!  (n  -  a;i  -  x2)!  (*2  -  ^2)!  (n  -  h  - t2 +xt  +  x2)\  t3\  (2n  -tl-  t3)\  ' 

Using  the  fact  that  xx  +  x2  +  x3  =  n  and  t1+t2+t3  =  2n,  this  ratio  reduces  to  1. 
The  result  clearly  holds  when  interchanging  the  other  columns  as  well. 
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3.4    The  Asymptotic  Approach  to  the  2  x  3  Contingency  Table 
The  asymptotic  approach  uses  the  Pearson  chi-square  statistic,  given  by 


2    3  /  _e..\2 

_  ^  ^  \ui]  C»J/ 
i=l  j=l  eij 


where  is  the  observed  value  in  the  ij  cell  and  is  the  expected  value  in  the  ij 
cell.  For  the  2  by  3  table 


■th 


x2 

X3 

yi 

V2 

V3 

with  equal  row  totals,  this  statistic  reduces  to 


{%i  Di) 


This  test  statistic  tends  to  be  liberal,  especially  for  small  sample  sizes,  as  was  seen  in 
the  examples  from  chapter  1. 

In  order  to  determine  the  necessary  sample  size  required  to  achieve  a  given  power, 
the  method  of  Lachin  (1977)  will  be  used.  This  method  is  based  on  the  non-central 
chi-square  disribution  and  is  described  as  follows.  For  all  chi-square  tests,  the  corre- 
sponding parameter  of  non-centrality,  call  it  A,  can  be  expressed  as  A  =  Nt(j1,'y2), 
where  N  is  the  total  sample  size  and  r  is  a  function  of  the  parameters  j1  and  j2 
which  apply  under  the  alternative  hypothesis.  For  a  chi-square  distribution  with  u 
degrees  of  freedom,  the  value  A(u,  a,  /?)  which  yields  power  (3  with  significance  level 
a  has  been  tabulated  by  Pearson  and  Hartley  (1972).  Given  the  required  A  and  the 
parameters  7X  and  72,  the  required  sample  size  is  obtained  as 


N  = 


r(7i,72)' 
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Thus,  all  that  is  needed  is  r(71,72).  For  the  2  by  3  table,  which  consists  of  two 
independent  trinomial  populations,  the  value  of  r(71,72)  is  given  by 


4  U  <*i 


where  71  j  and  j2j  are  the  values  of  the  parameters  for  populations  1  and  2  respectively 
under  the  alternative  hypothesis,  and 


1  2 


—  9     T»j  ■ 
L  »=i 


3.5   An  Exact  Unconditional  Test  for  the  2x3  Contingency  Table 

In  this  section,  the  methods  that  were  developed  in  chapter  2  are  used  to  determine 
the  size  of  the  Pearson  chi-square  test  statistic.  This  statistic  is  given  by 

2     3    /      _  \2 
j-<  _  ^2  ^  \uv  C'J/ 

i-l  j=l  eij 

where  Oy  is  the  observed  value  in  the  ijth  cell  and  is  the  expected  value  in  the  ijth 
cell.  The  asymptotic  null  distribution  of  T,  the  chi-square  distribution  with  2  degrees 
of  freedom,  is  typically  used  to  approximate  the  actual  size.  Hence,  the  results  of  the 
current  section  can  be  used  to  determine  the  accuracy  of  the  approximation. 

Since  X  and  Y  are  independent  trinomial  random  variables  with  parameters 
(n,  0X)  and  (n,  0y)  respectively,  the  power  function  is  given  by 


=  y(      U      )(      n      \f\ 6x,+y' 
Z£\x1x2x3)\y1  y2  y3  J  1 


44 


where  R  is  the  rejection  region  and  0,  the  nuisance  parameter,  is  the  common  value  of 
6X  and  6y  under  H0.  For  the  hypothesis  of  interest,  H0  :  0X  -  &y  versus  Ha  :  6X  ^  9y, 
the  form  of  the  critical  region  is  given  by 

R  =  {(x,y)  :  T  >  t,  xhyi  =  0,  ...,n,  £  x<  =  £  &  =  n}. 

For  the  purposes  of  exact  unconditional  inference,  we  need  to  find  a  value  of  T, 
call  it  Ta,  such  that 

Ta  =  inf{f  :  sup  71(61,62)  <  a} 
#ee0 

where 

60  =  {(0i,02)  :  0  <  0!,02  <  1,0  <  0i+02  <  1}. 

Before  using  the  methodology  of  chapter  2  to  determine  what  the  values  of  Ta  are, 
we  note  the  following  three  simplifications  that  can  be  made  to  reduce  the  amount  of 
calculations. 

1.  The  null  power  functions  corresponding  to  the  tables 


and 


Xi 

%2 

X3 

yi 

2/2 

2/3 

yi 

2/2 

2/3 

X2 

£3 

are  identical,  namely 

v{6u62)  =  y  (      n      )(     n     )f]6Xi+y'  . 
r\xix*x*  )  V  ^1  2/2  2/3  )l\  1 

Thus,  when  calculating  the  quantity  sup7r(0i,  02),  it  is  sufficient  to  consider  only  those 
tables  in  the  critical  region  which  are  not  row  permutations  of  any  other  table,  and 
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then  multiply  the  resulting  probability  by  2.  That  is,  the  null  power  function  can  be 
expressed  as 

where 

R*  =  {(x,y)  €  Vn  :  T  >  t,   xu  y{  =  0, n  ,  =  = 

and 

K  =  the  set  of  all  possible  tables  for  a  given  value  of  n,  such  that  any  table  belonging 
to  Vn  is  not  a  row  permutation  of  any  other  table  in  Vn. 
2.  The  nuisance  parameter  space,  given  by 

60  =  {(0i,02)  :  0  <  0i,02  <  1,0  <  0x  +02  <  1} 

is  the  set  over  which  the  size  of  the  test  is  to  be  determined.  This  set  can  be  re- 
duced, and  hence  make  calculations  easier,  by  noting  that  the  null  power  function  is 
symmetric  about  the  line  0X  —  02.  To  see  this,  consider  any  point  6  G  60,  such  that 
0i  >  02  (this  is  arbitrary  as  we  could  have  chosen  a  point  such  that  9\  <  02).  The 
null  power  function  at  such  a  point  is  given  by 


It  is  important  to  note  that  for  any  critical  region  of  the  form  R* , 

{{mxuyi,xi,yi):(x,y)eir}   =   {(mX2iW,  x2,  y2)  ■  (x,  y)  G  R*} 

=   {(jnX3,ya,X3,y3)  :  (x,y)  G  R*}  (3.1) 

where  mXuyi  is  the  number  of  such  (£i,?/i)  values.  This  follows  since,  for  a  given 
2  by  3  table,  the  Pearson  chi-square  test  statistic  is  invariant  to  row  and  column 
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permutations.  Thus, 


r*  V 


n     |  (     n     ]  ex,1+m6¥+V2(\  -  e,  -  e2)2n-xi-x?-yi-y* 

xi  x2x3  I  \  yi  y2  y3 


is  equal  to 


ar2  ar3  M  fi  ^  y3 


and  so  there  is  symmetry  about  the  line  9\  =  92.  Therefore,  we  need  only  consider 
the  null  power  function  over  the  region 

e0i^{(e1,92)eeQ:91>92} 

in  order  to  determine  the  size  of  the  test. 

3.  A  further  reduction  of  the  null  parameter  space  can  be  made  by  restricting  9\.  It 
will  be  shown  that  the  null  power  function  need  only  be  considered  over  the  set  60i 
with  the  additional  restriction  that  6\  <  0.5,  i.e.  the  set 

eO2  =  {(0i,02)Geoi:0i<o.5}  • 

This  again  follows  from  the  structure  of  the  critical  region  that  is  formed  based  on 

Pearson's  chi-square  statistic.  Consider  an  arbitrary  point,  (0i,02)  £  ©oi,  but  not  in 

002-  Define 

7!  =  1  -  9i  -  02,  and 

72  =  02- 

Then,  ji  <  0.5,  since  9X  +  02  >  0.5,  so  that  (71,72)  €  ©02-  Thus,  for  any  point 
(0i,  02)  £  ©oi /©02,  there  is  a  corresponding  point  (71, 72)  G  902  with  jx  =  1  -  Qx  -  02 
and  72  =  02 .  Then,  from  (3.1)  it  follows  that 


7^,02)  =  2  W      n      )(      n  )f[9?+« 
j£  \  xx  x2  x3  J  \  yx  y2  y3  J  11  « 
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is  equal  to 

*<«»=2?(I1:IJ)(,;Jr"' 

Thus,  for  any  point  in  601/602  there  exists  a  point  in  0O2  that  has  the  same  value  of 
the  null  power  function,  thus  we  need  only  consider  calculating  the  supremum  over 
the  region  0O2. 

In  summary,  the  problem  of  determining  the  size  of  the  Pearson  test  statistic 
amounts  to  maximizing  the  function 


over  the  set 

eoa  =  {(01,fc)€eOi:0i<o.5}. 

Once  the  values  of  TQ  have  been  determined  for  sufficiently  large  n,  the  power  can 
be  computed  for  a  given  size  a  and  values  of  Gx  and  9y  in  the  alternative  parameter 
space.  The  minimum  sample  size  required  to  attain  a  power  of  0  and  significance 
level  a  can  be  computed  by  solving 

ne  =  min{n  :  n(0x,  0y)  >  0} 

where  the  critical  region  that  defines  n(0x,  9y)  is  based  on  Ta,  a  function  of  n. 

3.6  Results 

For  n  =  10(1)70  and  a  =  0.025  and  0.05,  the  exact  critical  values,  TQ,  and  the 
size  (of  precision  0.001)  of  Ta  were  computed  and  given  in  Table  B.l  of  Appendix  B. 
Using  the  values  of  Ta  from  Table  B.l,  it  is  possible  to  compute  the  power  of  a  test 
of  level  a  =  0.025  or  0.05  for  given  values  of  0X  and  0y.  Sample  sizes  necessary  to 
attain  a  power  of  0  can  be  determined  by  solving 

ne  =  min{n  :  7r(0x,By)  >  0} 
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where  the  critical  region  is  determined  by  Ta,  which  depends  on  n.  This  equation 
was  solved  for  the  a  value  0.025  and  power  of  at  least  0.80.  The  results  are  given  in 
Table  B.2  of  Appendix  B  for  various  selected  values  of  0X  and  0y. 

A  methodology  was  developed  in  section  3.3  of  this  chapter  to  determine  the 
critical  regions  associated  with  FET,  which  in  turn  can  be  used  to  determine  the 
power  of  the  test,  and  hence  sample  size  requirements.  Table  B.3  in  Appendix  B 
gives  the  sample  sizes  required  to  achieve  a  power  of  at  least  0.80,  at  level  a  =  0.025, 
for  the  0X  and  0y  values  that  were  used  in  Table  B.2.  Thus,  comparisons  to  the  exact 
unconditional  method  can  be  made.  It  can  be  seen  that  neither  test  is  uniformly  best. 
In  many  cases  the  required  sample  sizes  are  the  same,  with  FET  statistic  achieving 
greater  power  in  the  majority  of  cases.  There  are  also  cases  were  FET  has  a  smaller 
required  sample  size,  and  one  case  of  the  exact  unconditional  method  requiring  a 
smaller  sample  size.  This  prompted  an  investigation  of  the  critical  regions  for  these 
two  testing  procedures.  For  a  —  0.025  and  n  =  10(1)40,  Table  B.4  in  Appendix  B 
lists  the  number  of  sample  points  in  the  critical  region  that  are  common  to  both  test 
statistics  and  the  number  of  points  that  are  unique  to  each.  For  small  values  of  n  it 
can  be  seen  that  the  unconditional  test  statistic  has  the  majority  of  unique  sample 
points  in  its  critical  region.  However,  as  n  gets  larger  (>  17),  it  is  FET  which  has  the 
larger  percentage  of  unique  points  in  the  critical  region.  It  seems  somewhat  counter 
intuitive  that  these  two  statistics  generate  such  different  critical  regions.  Therefore 
it  seems  appropriate  to  give  an  example  of  this  result.  The  following  example  is  for 
the  case  n  =  10.  The  table 


7 

0 

3 

6 

4 

0 

has  a  conditional  p-value  of  0.0217,  and  therefore  will  be  a  part  of  the  critical  region 
for  FET  at  level  0.025.  However, 
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which  is  not  part  of  the  critical  region  for  the  exact  unconditional  test  statistic  at  a 
level  0.025.  This  follows  from  the  fact  that  for  Ta  =  7.20,  the  null  power  function 
evaluated  at  the  point  9X  =  0.333  and  02  =  0.333  is  greater  than  0.025.  The  above 
table  belongs  to  the  critical  region  of  FET,  but  not  to  the  critical  region  of  the  exact 
unconditional  test  based  on  Pearsons  chi-square  statistic. 
Consider  the  table 


4 

1 

5 

5 

5 

0 

for  which 

2     3    /  0.._e..)2       I       42  55 

T  =  YY  ^  =  -  +  —  +  —  =  7.778 

Hrl       en  9     6  5 

which  belongs  to  the  critical  region  for  the  exact  unconditional  test  statistic  at  a 
level  0.025,  since  Ta  =  7.27  for  n  =  10.  However,  the  conditional  p-value  for  this 
table  is  0.0283,  and  so  it  will  not  be  a  part  of  the  FET  critical  region.  This  example 
demonstrates  that  there  are  indeed  dissimilarities  between  these  two  test  statistics. 
Table  B.5  in  Appendix  B  compares  the  results  of  Tables  B.2  and  B.3  to  the 

asymptotic  sample  sizes  found  by  using  the  method  of  Lachin.  Both  FET  and  the 
exact  unconditional  methods  required  sample  sizes  that  were  uniformly  smaller  than 
for  the  asymptotic  method,  for  the  values  of  0X  and  0y  that  were  considered. 

Appendix  C  contains  graphs  of  the  null  power  function  for  FET  and  the  exact 
unconditional  test  for  n  =  30, 40,  and  50,  and  a  =  0.025.  We  can  see  that  FET  is 
doing  quite  well.  That  is,  the  null  power  function  is  relatively  flat  and  attains  values 
near  a  for  most  of  the  null  parameter  space.  The  exact  unconditional  test  is  not  as 
well  behaved.  Its  graph  is  not  as  stable,  i.e.  flat,  and  does  not  achieve  values  near  a 
as  consistently  as  FET. 

In  conclusion,  the  results  of  this  chapter  demonstrate  that  it  is  not  the  case  that 
FET  is  uniformly  better  than  the  exact  unconditional  method  based  on  the  Pearson 
chi  square  test  statistic.  Nor  is  it  the  case  that  the  exact  unconditional  test  based  on 
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Pearsons  chi  square  test  is  uniformly  better  than  FET.  This  does  not  conflict  with 
the  result  at  the  end  of  chapter  1,  since  two  different  test  statistics  were  used.  The 
sample  size  requirements  for  the  exact  unconditional  test  and  FET  were  uniformly 
smaller,  for  the  examples  considered  here,  than  those  for  the  asymptotic  test  using 
the  method  of  Lachin.  Also,  the  graphs  of  the  null  power  functions  suggest  that  the 
conditional  p-value  should  be  used  instead  of  the  Pearson  chi-square  statistic  for  the 
exact  unconditional  method.  This  will  be  done  in  the  near  future. 


CHAPTER  4 
THE  2x2x2  CONTINGENCY  TABLE 


4.1  Introduction 


In  practice,  studies  can  be  designed  to  assess  the  association  between  two  random 
variables  X  and  Y  by  arranging  the  data  into  several  2x2  tables.  It  could  be  the 
case  that  the  possible  association  between  X  and  Y  is  of  considerable  interest,  and 
hence  there  are  several  investigators,  acting  independently  of  one  another,  studying 
the  association.  Also,  the  association  may  be  investigated  in  different  types  of  popu- 
lations, for  example  males  and  females.  Finally,  a  study  may  require  stratifying  the 
samples  being  compared  with  respect  to  variables  that  are  known  to  be  associated 
with  X  and  Y,  with  each  stratum  having  its  own  2x2  table. 

In  this  chapter  we  consider  the  case  of  two  2x2  tables.  Each  table  will  represent 
the  comparison  of  independent  binomial  random  variables.  It  is  also  assumed  that 
the  two  tables  are  independent  of  each  other.  For  example,  we  could  be  considering 
an  experiment  in  which  a  new  treatment  is  to  be  tested  for  its  efficacy  relative  to  the 
standard  treatment  being  used,  and  this  is  done  at  two  different  hospitals.  That  is, 
we  have  two  stations  where  clinical  trials  are  performed  to  compare  a  new  treatment 
versus  a  standard  treatment.  Because  our  main  interest  is  in  study  design,  where 
symmetry  seems  to  be  desired,  the  number  of  trials  for  each  of  the  4  binomial  distri- 
butions will  be  the  same.  The  data  for  such  an  experiment  can  be  represented  in  a 
2x2x2  table  as  follows: 


n  —  x\ 

y\ 

n-yi 

X2 

n  -  x2 

V2 

n-y2 
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52 


where  Xi,X2,Yi,  and  Y2  are  mutually  independent,  and 

X1  ~  Bi(n,6i) 

X2  ~  02 ) 
Yi  ~  Si(n,Ai) 

Y2  ~  fli(n,Aa). 
The  hypothesis  of  interest  is  formally  stated  as 


versus 


The  probability  of  observing  the  above  2x2x2  table  is 


where  jt  is  the  common  value  of  Q{  and  \  under  H0,  for  %  —  1, 2  . 

Because  of  the  nuisance  parameters  71  and  72,  asymptotic  or  conditional  methods 
have  traditionally  been  used  to  perform  the  above  hypothesis  test.  There  have  been  no 
attempts  to  compute  the  exact  unconditional  test  size,  because  of  the  computational 
complexity  involved,  and  hence,  the  maximization  method  of  Barnard(1947)  has  not 
been  used  for  this  hypothesis  testing  problem.  In  this  chapter,  methodology  developed 
in  chapter  2  will  be  used  to  compute  the  exact  unconditional  test  sizes,  allowing  us 
to  form  critical  regions  using  the  maximization  method.  First,  we  briefly  review  the 
conditional  and  asymptotic  approaches  to  this  problem. 
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4.2    The  Conditional  Approach  to  the  2  x  2  x  2  Contingency  Table 

The  joint  binomial  probability  distribution  is  an  exponential  family  distribution. 
Thus,  nuisance  parameters  can  be  eliminated  by  conditioning  on  appropriate  sufficient 
statistics.  After  conditioning,  the  resulting  probability  distribution  will  be  completely 
specified  under  H0.  The  joint  distribution  of  X\,X2,  Y\,  and  Y2  is  given  by 


^> =(:)(:)( : )  ( ; )  ft™  -*>-*  - 


D(x,  y)C(0,  A)  exp  j  J^fo  +  y{)  log  +  J^x<  log  '     _  ' 


where 


£>(x,y) 


nwnwnwn 
^i  /  V  ^2  /  V  2/i  /  V  2/2 


and 


c(e,\)  =  f[(i-9irii(i-\jr. 

i=i  j=\ 

The  above  probability  distribution  will  be  completely  specified  under  H0  if  we  con- 
dition on  Xi  +  Y\  and  X2  +  Y2.  Before  deriving  this  conditional  distribution,  define 
the  following 

t  =  {ti,t2)  i  where  ti  —     +  yi  for  i  =  1,  2 

At  =  [(x,  y) :     +  j/i  =  ti  i  =  1, 2  0  <  x,-,  ^  <  n  ] 

ai  =  max(0,  ti  —  n) 

a2  =  min(n,  ti) 

bi  =  max(0,  t2  —  n) 

b2  —  min(n,  t2). 

It  follows  that, 
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P(x|xi+yi  =  tux2  +  y2  =  h) 


D(x,  y)C(fl,  A)  exp  {E?=i  *i      (A)  +  Hi  ^  log 
D(x>,  y')C(0,  A)  exp  {eL  U  log  ^  +  E?=i  A  log  gg^j} 


^(x,y)exp{EL^log^}} 

EA^yO«p{5ii«5iogJ^Jj}} 


Finally,  since 


Il'=Ql  I2'=6l 


we  have  that 


P(X  =  x\X1+Yi  =  tuX2  +  Y2 
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where 

max(0,  h  —  n)  <  X\  <  min(n,  tx)  and, 
max(0,  t2  -  n)  <  x2  <  min(n,  t2). 

We  can  now  determine  the  critical  region  for  a  test  of  level  a.  Let  Po(x|t)  denote 
P(X  =  x\Xi  +  Yi  =  U,  X2  +  Y2  =  t2)  when  H0  is  true,  and  define 

p0(x*|t)  =  £P0(x|t) 

where 

Bt  =  [(x,y):P0(x|t)<P0(x*|t)]  . 

For  each  value  of  t,  let 

Cn,t(a)  =  [(x,y):po(x|t)<a)]  . 
Then,  from  Freeman  and  Halton  (1951),  the  critical  region  is  the  set 

Pn(a)  =  UC'«,t(«)  • 
t 

If  the  size  of  the  test  is  less  than  a,  which  will  usually  be  the  case  due  to  the  discrete 
null  distribution,  a  randomization  can  be  used  to  achieve  size  a.  The  randomized 
test  is  given  by 


0(x,t) 


1     x  G  Cnit(o!) 
7t   x  =  xt 
0  o/w 


where  xt  e  At\Ct(a)  and  p0(x(|t)  <  p0(x|t)  Vx  e  A\C„)t(a). 
Also, 

a-  Ecnt(Q)Po(x|t) 
7t  Po(x(|t) 

However,  randomized  tests  are  not  typically  used,  so  that  the  resulting  nonrandomized 
test  may  be  conservative. 
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4.3    Determining  Sample  Size  Requirements  Using  the  Conditional  Approach 


In  order  to  determine  the  sample  size  required  to  attain  a  given  power,  critical 
regions  must  be  determined  for  sufficiently  large  n.  As  was  the  case  for  the  2  x  3  table, 
the  critical  region  gets  quite  large  with  increasing  n.  This  causes  the  computing  time 
needed  to  form  critical  regions  for  even  moderate  n  to  be  much  too  great.  To  see  this, 
note  that  for  a  given  value  of  n,  the  possible  number  of  2  x  2  x  2  contingency  tables 
is 

[n  +  l]4  . 

A  simplification,  analogous  to  that  of  the  previous  chapter,  can  be  made  which 
reduces  the  calculations  and  amount  of  computer  time  considerably.  It  will  be  shown 
that  it  is  possible  to  calculate  the  power  using  a  condensed  sample  space.  This 
simplifies  calculations  and  reduces  the  amount  of  computer  time  needed  substantially. 
Define  the  following: 

Wn  =  the  set  of  all  2  x  2  tables  such  that  any  table  in  Wn  is  not  a  row  or  column 
permutation  of  any  other  table  in  Wn,  and  the  column  1  total  is  less  than  or  equal 

to  the  column  2  total. 
It  is  easy  to  show  that  for  a  given  value  of  n,  Wn  consists  of  the  tables 


i 

n  —  i 

3 

n  -  j 

where 


0  <  i  < 


i  <  j  <  n  —  i 


and  [*]  is  the  greatest  integer  function.  Sample  size  determinations  can  be  made  using 
only  those  2x2  tables  which  are  in  Wn,  which  will  be  demonstrated  next.  However, 
first  note  that  the  total  number  of  2  x  2  x  2  tables  based  on  W„  is 


L  2  J  n—i 

2 

[a] 

2 

EE 

X>-a+i) 

i=0  j=i 

i=0 

+  l)(n  +  l- 


n 
L2 
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which  equals 


(f +  1)4 


i2 


if  n  is  even 
if  n  is  odd 


and  note  that 


11  < 


(n  +  lY 


([2] +!)(„  +  !  _[*])] 


2  — 


<  16 


for  n  >  10  .  Hence,  for  reasonable  values  of  n,  the  sample  space  is  condensed  by  a 
factor  of  at  least  11. 

For  the  purposes  of  calculating  the  conditional  p- value  for  a  given  table,  and  hence 
determining  whether  or  not  it  will  be  in  the  critical  region,  we  must  consider  the  fol- 
lowing. When  calculating  the  conditional  p-value  for  a  given  table  considering  the 
entire  set  of  tables,  not  just  those  belonging  to  Wn,  the  usual  method  would  be  to 
first  calculate  the  conditional  probability  of  the  given  table, 


n 
•n 


n 


n 


n 


2n 


2n 


which  we  denoted  as  P0(x\t).  Then,  the  conditional  p-value  is  given  by 


P0((x/,y/)  6  At  :  P0(x/|t)  <  Po(x|t)). 


If  this  p-value  is  less  then  a  specified  a,  the  point  (x,  y)  will  be  included  in  the  critical 
region.  However,  when  using  only  those  tables  in  Wn,  the  probability  Po(x|t)  must 
be  modified.  The  modified  probabilities  will  be  denoted  by  P0*(x|t),  and  determined 
by  considering  the  possible  "types"  of  2  x  2  x  2  tables  that  can  be  generated  when 
X\,  X2,  Yi,  and  Y2  have  the  same  number  of  trials,  n.  Define  the  following  types  of 
2x2  tables: 
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unique:  A  table  is  unique  if  no  cells  are  the  same,  for  example 


0 

5 

3 

2 

diagonal  symmetric:  A  table  is  diagonal  symmetric  if  Xi  =  n  -  yi  ^  n  -  X(,  for 
example 


0 

5 

5 

0 

symmetric:  A  table  is  symmetric  if  Xi  =  yt,  for  example 


Now,  the  2x2x2  table  types  and  associated  probability  values  P0*(x|t)  are  defined 
as  follows: 

a.  Both  tables  are  unique.  For  example, 


0 

5 

3 

2 

1 

4 

2 

3 

In  this  case  P0*(x|t)  =  4  x  P0(x|t). 

b.  One  table  is  unique  and  the  other  is  symmetric.  For  example, 


0 

5 

4 

1 

In  this  case  P0*(x|t)  =  2  x  P0(x|t). 

c.  One  table  is  unique  and  the  other  is  diagonal  symmetric.  For  example, 
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0 

5 

4 

1 

0 

5 

5 

0 

In  this  case  P0*(x|t)  =  4  x  P0(x|t). 

d.  Both  tables  are  symmetric.  For  example, 


In  this  case  P0*(x|t)  =  Po(x|t). 

e.  One  table  is  symmetric  and  the  other  is  diagonal  symmetric.  For  example, 


In  this  case  P0*(x|t)  =  2  x  P0(x|t). 

f.  Both  tables  are  diagonal  symmetric.  For  example, 


0 

4 

4 

0 

0 

4 

4 

0 

In  this  case  P0*(x|t)  =  4  x  P0(x|t). 
This  exhausts  the  possibilities. 

The  p- value  for  any  table  with  entries  x\,  X2,  y\,  and  y^  is  calculated  as 


where 


P((x/,y/)e^P*(x/|t)<P0*(x|t)) 


A*t  =  (Wn®Wn)f)At. 


The  critical  region  will  contain  all  those  tables  in  Wn  (g)  Wn  that  have  corresponding 
p-values  <  a.  Once  the  critical  regions  are  formed  for  a  given  a  level,  it  is  possible 
to  calculate  the  power.  First,  define  the  following 


Ci(0i,02,  Ai,  a2)  =  e\e{\kxxs2(\  -  00^(1  -  02)n-J'(i  -  A!)n-*(i  -  a2)"- 
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51  =  {{xi,yi,x2, 1/2),  (xi,yi,y2, x2),  {xuVun-  y2,  n  -  x2),  (xi, yu n  -  x2,  n  -  y2), 
{yi,xux2,  y2),  (yi,xu  y2,  x2),  (yux^n  -y2,n-  x2),  {y\,xun-  x2,  n-y2), 

(n -xi,n- yi,x2,y2),  (n  -  xun  -  yi,y2,x2),  {n  -  xun  -  yu  n  -  y2,n-  x2), 
(n  -Xun-  yi,n-  x2,n-  y2),  (n  -yun-  xux2,y2),  (n  -yun-  xuy2,x2), 
(n  -yi,n-  xun-  y2,n-  x2),  (n  -yun-  xun-  x2,n-  y2)} 

52  =  {(xi,yi,x2,y2),  (xi,yi,y2,x2),  (yi,xux2,y2),  {yi,xx,y2,x2), 

(n  -xi,n-  yi,x2,y2),  (n  -xun-  yi,y2,x2),  (n  -yun-  xux2,y2), 
(n  -yx,n-  xx,y2,x2)} 

53  =  {{xi,yi,x2,y2),  (xi,yuy2,x2),  (yi,xux2,y2),  (yi,xi,y2,x2), 
(xuyun-  x2,n-  y2),  (xuyun-  y2,n-  x2),  {yuxun-  x2,n-  y2), 
(yuxun-  y2,n-  x2)} 

Sa  =  {(xi,yi,x2,y2),(xuyi,n-  x2,n-  y2),  (yu  xu  x2,  y2),  (yu  xu  n  -  x2,  n  -  y2), 
(n  -xi,n-  yi,x2,y2),  (n  -xun-  yun-  x2,n-  y2),  (n  -yun-  xx,x2,y2), 
(n  -y\,n-  xun-  x2,n-  y2)} 

S5  =  {(xuVi,x2,  y2),  (xuyi,  y2,  x2),  {xuyun-  x2,  n  -  y2),  (xu  yu  n  -  y2,  n-x2), 
(n  -xun-  yi,x2,y2),  (n -xx,n- yuy2,x2),  (n  -  xun  -  yun-  x2,n-  y2), 
(n  -xun-  yun-  y2,n-  x2)} 

Se  =  {(xuyux2,y2),  (xi,yuy2,x2),  (yuxi,x2,y2),  (yi,xuy2,x2)} 

S7  =  {(xi,yi,x2,y2),  (n  -  xun  -  yux2,y2),  (yi,xux2,y2),  (n  -yun-  xux2,y2)} 

Ss  =  {(xi,yi,x2,y2),  (xuy1,y2,x2),  (xuyun-  x2,n-  y2),  (xx,yx,n-  y2,n-  x2)} 

S9  =  {(xi,yi,x2,y2),  (n  -  xun  -  yux2,y2),  (xuyi,n-  x2,n-  y2), 
(n  -xi,n-  yun-  x2,n-  y2)} 


>10 


5ii  = 


>5l2  — 


s 
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5l5  = 


«i,  l/i,  X2, 2/2),  (j/i,  xi.flfa,  2/2),  yi,n-x2,n-  y2),  (j/i,  *i,  n  -  x2,  n  - 
xi,yi,x2,  jfe),  (xi,  1/1, 2/2,  z2),  (n  -  xi,  n  -  yux2,  y2),(n-  x^n-  yuy2, 


xi,  yi,x2,y2),  {xuyi,n-  x2,n-  y2)} 


xi,yi,x2,y2),  (n  -  xun  -  yi,x2,y2)} 


xi,yi,x2,y2),(xuyi,y2,x2)} 


x\,yi,x2,y2),  (yi,xux2,y2)} 


Sw  =  {{x\,yi,x2,y2)} . 


The  power  function  is  comprised  of  the  following  types  of  terms, 

1.  For  contingency  tables  of  the  form  (a)  the  contribution  to  the  power  is 

D(x,y)£Ci(0i,fc,Ai,A2). 

Si 

2.  For  contingency  tables  of  the  form  (b)  the  contribution  to  the  power  is 

D(x,y)£C1(0i,fclAi,Aa) 
s4 


or 


or 


or 


D(x,y)  02,  Al5A2) 


D(x,y)£Ci(0i,fc,Ai,A2) 
s7 


D(x,y)j;C1(tf1>fl2,Ai,A2). 

s8 


3.  For  contingency  tables  of  the  form  (c)  the  contribution  to  the  power  is 


£>(x,y)  ^d^^a.A^Aa) 
s2 
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or 


JD(xfy);£C7i(0i,ft)Ai,Aa). 

4.  For  contingency  tables  of  the  form  (d)  the  contribution  to  the  power  is 

D(x,y)£Ci(*i,fc,Ai.Aa) 
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or 


or 


or 


D(x,y)Y^cl(el,e2,xux2) 

Sl6 


D(x,y)Ylcl(eue2,Xu\2) 

Sl2 


s13 


5.  For  contingency  tables  of  the  form  (e)  the  contribution  to  the  power  is 

D(x,y)£Ci(tf1,ft,A1>AI) 

Sit 

or 

D(x,y)X;Ci(^i^2)Ai,A2). 

6.  For  contingency  tables  of  the  form  (f)  the  contribution  to  the  power  is 

D(x,y)^C1(^1^2,A1,A2). 

Now,  proceed  by  iterations  on  n  to  determine  the  sample  size  required  to  achieve  the 
specified  power.  This  can  be  expressed  as 

nc  =  inf{n  :  nn(0u62,  Xu  X2)  >  /?} 

where  the  subscript  c  denotes  the  sample  size  is  based  on  conditional  methods,  @  is 
the  desired  power,  and  nn(6u  92,  Xu  A2)  represents  the  power  function  for  a  given  n. 
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4.4    The  Asymptotic  Approach  to  the  2  x  2  x  2  Contingency  Table 

In  what  follows,  it  is  convenient  to  think  of  each  2x2  table  as  representing  the 
comparison  of  a  study  and  control  group.  Define 
N  :  the  total  number  of  observations 

riijk  :  the  number  of  observations  in  the  jth  group  of  the  ith  table  with  outcome  k. 
Vi  :  proportion  of  observations  in  table  i. 

Pi  :  proportion  of  observations  in  the  study  group  in  the  ith  table. 
iTij  :  probability  of  success  for  the  jth  group  in  the  ith  table. 

TTi  :  Pi^il  +  (1  -  Pi)^i2 
Si  :  ITn  -  7T,2 

Q     _  7Til(l-Ti2) 
»  ^(l-Xil) 

where  i,  j,  A;  =  1,2.  Also,  j  =  1  corresponds  to  the  study  group,  and  k  =  1  corresponds 
to  a  success. 

We  consider  two  forms  of  the  Mantel-Haenszel  statistic, 
1.  without  the  continuity  correction 

M2  =  (EM2 


2.  with  the  continuity  correction 


M. 


c  ^2 


where 

Pi  =  NuiPi(l  -  pi)(pn  -  Pa) 
y.  _  "ti  "i2  Pi(l-Pt) 

Pi?  =       and  pi  =  rh^L  ■ 

r  J         Uij.  rl  m.. 

These  statistics  have  asymptotic  chi-square  distributions  with  1  degree  of  freedom 
under  H0  :  9^  =  02  =  1 . 
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To  determine  sample  size  requirements,  results  from  Wittes  and  Wallenstein(1987) 
will  be  used.  When  the  number  of  tables,  T,  is  fixed,  and  nj..  — >•  oo  in  such  a  way 
that  y/NSi  approaches  a  nonzero  constant  7$,  they  have  demonstrated  the  following 

d,   jvr,..  2 


Nfao*)  (4.1) 


EV>    P+a2  (4.2) 


'N 

where  fx  and  a2  are  given  by 

i 

and 

nij  — >  Ki    as    nj..  — >  00. 

From  (4.1),  (4.2),  and  a  theorem  of  Slutsky,  the  power  of  the  Mantel-Haenszel 
statistic  without  the  continuity  correction  is  given  by 

p,  {m2  >  *?(«)}  =  p„  {  J^l  >  y^J + p,  {J^  <  -y^wj 
=  *{§-zt}  +  *{-|-zs}  +  <*(i), 

and  the  power  of  the  Mantel-Haenszel  statistic  with  continuity  correction  is  given  by 

p,  {mi  >  Xfw}  =  p„  JE*y  >        +  P<i  {Eg_d  < 

-.{-•z£-^}  +  .{-£+£_zf}  +  %(1, 


where 

1 


a  = 


2\/N 


and  a2  =  ^  ViPi{\  -  Pi)iXi{\  -  7^). 
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Therefore,  for  a  given  level  a,  the  sample  size  requirements  are  determined  for 
each  statistic  by  calculating  the  power  iteratively  with  respect  to  N  until  the  desired 
P  is  achieved. 

4.5    An  Exact  Unconditional  Test  for  the  2  x  2  x  2  Contingency  Table 

The  methods  that  were  developed  in  chapter  2  will  be  used  to  determine  the  size 
of  the  Mantel-Haenszel  test  statistic.  This  statistic  is  given  by 

2  _  I  (Zj=i9i\  ~  j) 

c  ~     £L  vt 

where  and  V*  were  given  in  the  previous  section.  The  asymptotic  null  distribution 
of  Ml ,  the  chi-square  distribution  with  1  degree  of  freedom,  is  typically  used  to  ap- 
proximate the  actual  size.  The  results  of  the  current  section  can  be  used  to  determine 
the  accuracy  of  the  approximation. 

Since  X\,  X2,  Yi,  and  Y2  are  independent  binomial  random  variables  with  pa- 
rameters (n,  6\),  (71,62),  (n,Xi),  and  (n,  A2)  respectively,  the  power  function  is  given 
by 

7r(0!,02,  Ax,A2)  =  ]>>(x, y)  f[0?*(l  -  6^  {[  Af  (1  -  A*)""" 

R  i=l  i=l 

=  E^y)n7r+y,(i-7,)2n-ii-yt 

R  i=l 

where  R  is  the  rejection  region,  is  the  common  value  of  6X  and  Ai  under  H0,  and  72 
is  the  common  value  of  62  and  A2  under  H0.  The  form  of  the  critical  region  is  given 
by 

R  =  {(x,  y)  :  M2C  >  m,  xu  y{  =  0, n). 

For  the  purposes  of  exact  unconditional  inference,  we  need  to  find  a  value  of  Mc2,  call 
it  Ma,  such  that 

Ma  =  inf{m  :  sup  ^(71,72)  <  a} 
7er0 
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where 


r„  =  {(7i,72)  :0<  7i,72  <!}• 


Before  using  the  methodology  of  chapter  2  to  determine  what  the  values  of  Ma 
are,  we  note  the  following  two  simplifications  that  can  be  made  to  reduce  the  amount 
of  calculations. 

1.  The  null  power  functions  corresponding  to  the  tables 


Xi 

n  —  X\ 

y\ 

n-yi 

X2 

n  -  x2 

V2 

n-yi 

and 


y\ 

n-yi 

n  —  x\ 

V2 

n-V2 

X2 

n  -  x2 

are  identical.  Thus,  when  calculating  the  quantity  sup 7r(7i,  72),  it  is  sufficient  to 
consider  only  those  2x2x2  tables  in  the  critical  region  which  remain  unique  upon 
the  permutation  of  both  rows.  That  is,  the  null  power  function  can  be  expressed  as 


tt(7i,  72)  =  2  £  D(x,  y)  II  7?i+w (1  "  702""*"* 


2 

[I 


where 


R* 


{(x,y)  G  Sn  :  Mc2  >  m,  xhyi  =  0, n},  and 


Sn  =  the  set  of  all  possible  2x2x2  tables  for  a  given  value  of  n,  such  that  any  table 
belonging  to  Sn  is  not  a  permutation,  wrt  to  rows,  of  any  other  table  in  Sn. 
2.  The  null  parameter  space,  given  by 

ro  =  {(7i,72):0<^,^2<l}, 

is  the  set  over  which  the  size  of  the  test  is  to  be  determined.  Using  the  form  of  the 
null  power  function,  and  the  fact  that  the  Mantel-Haenszel  statistic  is  invariant  to 
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row  and  column  permutations,  the  region  over  which  the  size  is  determined  can  be 
reduced  to 

Ti  =  {(71,72)  :  0  <  71  <  0.5  and72  <  7l}, 

thus, 

Ma  =  inf{ra  :  sup  ^(71,72)  <  a}. 

Once  the  values  of  Ma  have  been  determined  for  sufficiently  large  n,  the  power 
can  be  computed  for  a  given  level  a  and  values  of  9\,  92,  X\,  A2  in  the  alternative 
parameter  space.  The  minimum  sample  size  required  to  attain  a  power  of  (3  and 
significance  level  a  can  be  computed  by  solving 

ne  =  mm{n  :  77(81,62,  \i,  A2)  >  P} 

where  the  critical  region  that  defines  tx(9\,92,  \\,  A2)  is  based  on  Ma,  a  function  of  n. 

4.6  Results 

For  n  =  10(1)50  and  a  =  0.025  and  0.05,  the  exact  critical  values,  MQ,  and  the 
size  (of  precision  0.001)  of  Ma  were  computed  and  given  in  Table  D.l  of  Appendix  D. 
It  should  be  noted  that  the  exact  unconditional  critical  values,  Ma,  were  considerably 
larger  than  the  asymptotic  value,  even  for  large  values  of  n.  This  suggests  that  the 

approximation  Mc2  ~  x?  maY  n°t  be  very  good. 

Using  the  values  of  Ma  from  Table  D.l,  it  is  possible  to  compute  the  power  for 
a  —  0.025  or  0.05  and  given  values  of  #i,02,Ai,  and  A2.  Sample  sizes  necessary  to 
attain  a  power  of  (3  can  be  determined  by  solving 

ne  =  min{n  :  7r(0i ,  6>2,  Ai ,  A2)  >  0} 

where  the  critical  region  is  determined  by  Ma,  which  depends  on  n.  This  equation 
was  solved  for  a  =  0.025  and  a  power  of  at  least  0.80.  The  results  are  given  in  Table 
D.2  of  Appendix  D  for  various  selected  values  of  9\,  92,  \\,  and  A2. 
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A  methodology  was  developed  in  section  4.3  of  this  chapter  to  determine  the 
critical  regions  associated  with  FET,  which  in  turn  can  be  used  to  determine  the 
power,  and  hence  sample  size  requirements.  Table  D.3  in  Appendix  D  gives  the 
sample  sizes  for  the  61,62,  Xi,  and  A2  values  that  were  used  in  Table  D.2,  in  order 
to  make  comparisons  with  the  exact  unconditional  method.  It  can  be  seen  that 
neither  test  is  uniformly  best.  In  the  majority  of  cases  the  required  sample  sizes 
were  different  for  the  two  test  statistics,  with  FET  having  an  appreciable  advantage 
in  several  instances.  This  prompted  an  investigation  of  the  critical  regions  for  these 
two  testing  procedures.  For  a  =  0.025  and  n  —  10(1)50,  Table  D.4  in  Appendix  D 
lists  the  number  of  sample  points  in  the  critical  region  that  are  common  to  both  test 
statistics  and  the  number  of  points  that  are  unique  to  each.  It  can  be  seen  that  the 
FET  statistic  has  the  majority  of  unique  sample  points  in  its  critical  region  for  all 
but  one  value  (n  =  11). 

Appendix  E  contains  graphs  of  the  null  power  function  for  FET  and  the  exact 
unconditional  test  for  n  =  10,20,30,40,  and  50,  and  a  =  0.025.  We  can  see  that 
the  null  power  function  for  FET  is  considerably  more  stable  than  that  of  the  exact 
unconditional  method  based  on  the  Mantel  Haenszel  statistic.  That  is,  the  null  power 
function  is  flatter  and  attains  values  near  a  more  frequently.  The  exact  unconditional 
test  is  not  well  behaved.  Its  graph  is  not  stable,  as  there  appears  to  be  a  protuberance 
in  the  center.  This  explains  why  the  exact  unconditional  critical  values,  i.e.  the  MQ's, 
were  so  much  larger  than  the  asymptotic  values. 

In  conclusion,  the  results  of  this  chapter  demonstrate  that  it  is  not  the  case  that 
FET  is  uniformly  better  than  the  exact  unconditional  method  based  on  the  Mantel 
Haenszel  test  statistic.  Nor  is  it  the  case  that  the  exact  unconditional  test  based 
on  the  Mantel  Haenszel  test  statistic  is  uniformly  better  than  FET.  This  does  not 
conflict  with  the  result  at  the  end  of  chapter  1,  since  two  different  test  statistics 
were  used.  Also,  the  graphs  of  the  null  power  functions  suggest  that  the  conditional 
p-value  should  be  used  instead  of  the  Mantel  Haenszel  test  statistic  for  the  exact 
unconditional  method.  This  will  be  done  in  the  near  future. 


APPENDIX  A 
THE  SET  U3 


Here  are  the  elements  of  the  set  U3  : 


3 

0 

0 

3 

0 

0 

3 

0 

0 

3 

0 

0 

3 

0 

0 

3 

0 

0 

3 

0 

0 

2 

1 

0 

1 

2 

0 

0 
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0 

1 

1 

1 

0 

2 
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2 

1 
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2 

1 
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2 

1 

0 
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2 
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2 
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2 

1 
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1 

1 

1 

1 

2 

0 

1 

0 

1 

2 

1 

1 

1 

1 

1 

1 

with  3,  12,  12,  6,  6,  12,  6,  6,  6,  12,  12,  6,  and  1  representations,  respectively. 
For  example,  the  table 


2 

1 

0 

0 

1 

2 

represents  the  following  6  tables: 


2 

1 

0 

0 

1 

2 

1 

2 

0 

1 

0 

2 

2 

0 

1 

0 

2 

1 

0 

1 

2 

2 

1 

0 

1 

0 

2 

1 

2 

0 

0 

2 

1 

2 

0 

1 
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APPENDIX  B 
CHAPTER  3  TABLES 


The  following  are  legends  to  be  used  with  the  tables  in  this  appendix. 
Table  B.l 

n     Sample  size  for  each  trinomial  population. 

TQ  Boundary  of  the  critical  region  associated  with  the  exact  unconditional  test. 
a*   Approximate  (to  within  0.001)  size  of  the  test  based  on  Ta. 

Table  B.2 

ne    Sample  size  required  for  the  exact  unconditional  test  statistic. 

Ta   Boundary  of  the  critical  region  associated  with  the  exact  unconditional  test. 

P     Power  for  given  ne  and  parameter  values  6X  and  0y. 

Table  B.3 

nc    Sample  size  required  for  FET  statistic. 

/3     Power  for  given  nc  and  parameter  values  0X  and  6y. 

Table  B.4 

n  Sample  size  for  each  trinomial  population. 

Nn  Number  of  sample  points  for  given  n. 

Nt,  Number  of  critical  region  points  in  common. 

Nc  Number  of  points  unique  to  FET  critical  region. 

Ne  Number  of  points  unique  to  exact  unconditional  critical  region. 

Pe  The  proportion  NN*N  . 
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Table  B.5 

ne    Sample  size  required  for  the  exact  unconditional  test  statistic. 
na    Sample  size  required  for  the  asymptotic  test  statistic. 
nc    Sample  size  required  for  FET  statistic. 
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Table  B.l:  Ta  and  a*  values  for  the  comparison  of  two  independent  multinomial 
distributions. 


a  — 

0.025 

a  =  0.05 

n 

T 

a* 

n 

T 

a* 

10 

7.27 

0.02479 

10 

5.90 

0.04917 

11 

7.27 

0.02501 

11 

6.13 

0.04986 

12 

7.28 

0.02478 

12 

6.27 

0.04933 

13 

7.39 

0.02462 

13 

6.11 

0.04989 

14 

7.40 

0.02472 

14 

6.07 

0.04990 

15 

7.36 

0.02485 

15 

6.11 

0.04983 

16 

7.48 

0.02491 

16 

6.15 

0.04958 

17 

7.48 

0.02498 

17 

6.11 

0.04945 

18 

7.47 

0.02494 

18 

6.07 

0.04852 

19 

7.40 

0.02493 

19 

6.02 

0.04992 

20 

7.38 

0.02482 

20 

6.04 

0.04966 

21 

7.42 

0.02489 

21 

6.11 

0.04977 

22 

7.41 

0.02494 

22 

6.08 

0.04942 

23 

7.37 

0.02502 

23 

6.08 

0.04979 

24 

7.39 

0.02498 

24 

6.05 

0.04964 

25 

7.37 

0.02490 

25 

6.01 

0.04963 

26 

7.38 

0.02475 

26 

6.07 

0.04972 

27 

7.44 

0.02496 

27 

6.10 

0.04968 

28 

7.42 

0.02493 

28 

6.10 

0.04962 

29 

7.38 

0.02496 

29 

6.07 

0.04979 

30 

7.38 

0.02498 

30 

6.03 

0.05001 

31 

7.39 

0.02492 

31 

6.03 

0.04999 

32 

7.38 

0.02487 

32 

6.05 

0.04959 

33 

7.40 

0.02499 

33 

6.07 

0.04953 

34 

7.42 

0.02499 

34 

6.04 

0.04960 

35 

7.42 

0.02486 

35 

6.04 

0.04953 

36 

7.42 

0.02481 

36 

6.04 

0.04958 

37 

7.39 

0.02490 

37 

6.04 

0.04951 

38 

7.40 

0.02494 

38 

6.04 

0.04991 

Table  B.l  (con't) 


a  = 

0.025 

a  =  0.05 

n 

T 

T 

n 

T 

a* 

39 

7.43 

0.02488 

39 

6.05 

0.04989 

40 

7.45 

0.02492 

40 

6.08 

0.04967 

41 

7.43 

0.02492 

41 

6.05 

0.04966 

42 

7.44 

0.02491 

42 

6.06 

0.04950 

43 

7.42 

0.02485 

43 

6.06 

0.04990 

44 

7.43 

0.02486 

44 

6.06 

0.04955 

45 

7.46 

0.02498 

45 

6.04 

0.04995 

46 

7.44 

0.02499 

46 

6.05 

0.04982 

47 

7.43 

0.02489 

47 

6.05 

0.05000 

48 

7.44 

0.02493 

48 

6.05 

0.05001 

49 

7.44 

0.02501 

49 

6.04 

0.04982 

50 

7.44 

0.02499 

50 

6.04 

0.04996 

51 

7.43 

0.02491 

51 

6.05 

0.04969 

52 

7.43 

0.02499 

52 

6.04 

0.04982 

53 

7.44 

0.02487 

53 

6.04 

0.04995 

54 

7.43 

0.02499 

54 

6.03 

0.04976 

55 

7.43 

0.02496 

55 

6.04 

0.04994 

56 

7.44 

0.02500 

56 

6.06 

0.04971 

57 

7.44 

0.02492 

57 

6.04 

0.04992 

58 

7.44 

0.02491 

58 

6.04 

0.04984 

59 

7.43 

0.02495 

59 

6.04 

0.04966 

60 

7.44 

0.02498 

60 

6.05 

0.04949 

61 

7.44 

0.02492 

61 

6.05 

0.04962 

62 

7.44 

0.02496 

62 

6.04 

0.04960 

63 

7.43 

0.02493 

63 

6.03 

0.04999 

64 

7.43 

0.02501 

64 

6.04 

0.04993 

65 

7.43 

0.02499 

65 

6.05 

0.04967 

66 

7.44 

0.02494 

66 

6.04 

0.04982 

67 

7.43 

0.02491 

67 

6.03 

0.04977 

68 

7.43 

0.02492 

68 

6.03 

0.04984 

69 

7.44 

0.02495 

69 

6.03 

0.04997 

70 

7.43 

0.02492 

70 

6.03 

0.04984 
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Table  B.2:  Sample  sizes  necessary  to  achieve  a  power  of  at  least  0.80,  at  level 
a  =  0.025,  for  comparing  two  independent  multinomial  distributions  using  exact 
unconditional  methods. 


flu 

#12 

^21 

022 

ne 

T 

* 

a 

0 

0.05 

0.05 

0.50 

0.05 

20 

7.38 

0.02482 

0.8180 

0.05 

0.05 

0.35 

0.20 

23 

7.37 

0.02502 

0.8120 

0.05 

0.05 

0.05 

0.35 

36 

7.42 

0.02481 

0.8086 

0.05 

0.05 

0.20 

0.20 

46 

7.44 

0.02502 

0.8073 

0.05 

0.20 

0.05 

0.65 

25 

7.37 

0.02490 

0.8029 

0.05 

0.20 

0.35 

0.20 

35 

7.42 

0.02486 

0.8126 

0.05 

0.20 

0.05 

0.50 

55 

7.43 

0.02496 

0.8055 

0.05 

0.20 

0.20 

0.05 

61 

7.44 

0.02492 

0.8091 

0.05 

0.35 

0.65 

0.20 

12 

7.28 

0.02478 

0.8218 

0.05 

0.35 

0.20 

0.05 

32 

7.38 

0.02487 

0.8165 

0.05 

0.35 

0.05 

0.05 

36 

7.42 

0.02481 

0.8086 

0.05 

0.35 

0.05 

0.65 

59 

7.43 

0.02495 

0.8045 

0.20 

0.35 

0.65 

0.05 

22 

7.41 

0.02494 

0.8238 

0.20 

0.35 

0.50 

0.05 

30 

7.38 

0.02498 

0.8141 

0.20 

0.35 

0.50 

0.35 

40 

7.45 

0.02492 

0.8097 

0.20 

0.35 

0.50 

0.20 

56 

7.44 

0.02500 

0.8024 

0.35 

0.35 

0.05 

0.20 

23 

7.37 

0.02502 

0.8050 

0.35 

0.35 

0.05 

0.65 

33 

7.40 

0.02499 

0.8051 

0.35 

0.35 

0.05 

0.50 

37 

7.39 

0.02490 

0.8122 

0.35 

0.20 

0.65 

0.20 

46 

7.44 

0.02499 

0.8034 

0.35 

0.50 

0.50 

0.05 

19 

7.40 

0.02493 

0.8192 

0.35 

0.50 

0.05 

0.50 

28 

7.42 

0.02493 

0.8060 

0.35 

0.50 

0.05 

0.80 

35 

7.42 

0.02486 

0.8047 

0.35 

0.50 

0.65 

0.20 

51 

7.43 

0.02491 

0.8028 
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Table  B.3:  Sample  sizes  necessary  to  achieve  a  power  of  at  least  0.80,  at  level 
a  =  0.025,  for  comparing  two  independent  multinomial  distributions  using  exact 
unconditional  methods  using  conditional  methods. 


0ii 

#12 

^21 

022 

nc 

P 

0.05 

0.05 

0.50 

0.05 

19 

0.8277 

0.05 

0.05 

0.35 

0.20 

23 

0.8182 

0.05 

0.05 

0.05 

0.35 

34 

0.8078 

0.05 

0.05 

0.20 

0.20 

46 

0.8020 

0.05 

0.20 

0.05 

0.65 

23 

0.8006 

0.05 

0.20 

0.35 

0.20 

35 

0.8093 

0.05 

0.20 

0.05 

0.50 

52 

0.8004 

0.05 

0.20 

0.20 

0.05 

61 

0.8050 

0.05 

0.35 

0.65 

0.20 

13 

0.8388 

0.05 

0.35 

0.20 

0.05 

31 

0.8015 

0.05 

0.35 

0.05 

0.05 

34 

0.8078 

0.05 

0.35 

0.05 

0.65 

57 

0.8052 

0.20 

0.35 

0.65 

0.05 

22 

0.8134 

0.20 

0.35 

0.50 

0.05 

30 

0.8076 

0.20 

0.35 

0.50 

0.35 

40 

0.8070 

0.20 

0.35 

0.50 

0.20 

56 

0.8000 

0.35 

0.35 

0.05 

0.20 

24 

0.8181 

0.35 

0.35 

0.05 

0.65 

33 

0.8023 

0.35 

0.35 

0.05 

0.50 

37 

0.8040 

0.35 

0.20 

0.65 

0.20 

46 

0.8020 

0.35 

0.50 

0.50 

0.05 

19 

0.8048 

0.35 

0.50 

0.05 

0.50 

28 

0.8038 

0.35 

0.50 

0.05 

0.80 

36 

0.8153 

0.35 

0.50 

0.65 

0.20 

51 

0.8042 
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Table  B.4:  A  comparison  of  the  critical  regions  of  FET  and  the  unconditional 
method. 


n 

Nn 

Nb 

Nc 

Pe 

10 

4356 

1458 

18 

96 

84.2 

11 

6084 

2202 

60 

120 

66.7 

12 

8281 

3264 

60 

132 

68.8 

13 

11025 

4374 

60 

180 

75.0 

14 

14400 

6528 

186 

192 

50.8 

15 

18496 

8850 

138 

228 

62.3 

16 

23409 

11628 

300 

240 

44.4 

17 

29241 

15174 

258 

288 

52.7 

18 

36100 

19344 

432 

348 

44.6 

19 

44100 

24450 

522 

372 

41.6 

20 

53361 

30588 

558 

324 

36.7 

21 

64009 

37632 

576 

468 

44.8 

22 

76176 

45918 

630 

492 

43.9 

23 

90000 

55560 

600 

480 

44.4 

24 

105625 

66288 

918 

612 

40.0 

25 

123201 

78984 

810 

612 

43.0 

26 

142884 

93132 

1014 

708 

41.1 

27 

164836 

109272 

1068 

780 

42.2 

28 

189225 

127104 

1308 

828 

38.8 

29 

216225 

147300 

1302 

900 

40.9 

30 

146016 

169860 

1392 

936 

40.2 

31 

278784 

194616 

1644 

1260 

43.4 

32 

314721 

222300 

1746 

1560 

47.2 

33 

354025 

252762 

1758 

1704 

49.1 

34 

396900 

286104 

2076 

1584 

43.3 

35 

443556 

322638 

2280 

1620 

41.5 

36 

494209 

326544 

2622 

1836 

41.2 

37 

548081 

406578 

2358 

2076 

46.8 

38 

608400 

453978 

2352 

2388 

50.4 

39 

672400 

505128 

3114 

2472 

44.3 

40 

741321 

560922 

3240 

2412 

42.7 
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Table  B.5:  Comparisons  among  exact  unconditional,  asymptotic,  and  FET  sample 
sizes  for  comparing  two  independent  multinomial  distributions  at  level  a  —  0.025. 


On 

#21 

ne 

na 

nc 

0.05 

0.05 

0.50 

0.05 

20 

22 

19 

0.05 

0.05 

0.35 

0.20 

23 

25 

23 

0.05 

0.05 

0.05 

0.35 

36 

40 

34 

0.05 

0.05 

0.20 

0.20 

46 

48 

46 

0.05 

0.20 

0.05 

0.65 

25 

27 

23 

0.05 

0.20 

0.35 

0.20 

35 

38 

35 

0.05 

0.20 

0.05 

0.50 

55 

56 

52 

0.05 

0.20 

0.20 

0.05 

01 

64 

61 

0.05 

0.35 

0.65 

0.20 

12 

14 

13 

0.05 

0.35 

0.20 

0.05 

32 

35 

31 

0.05 

0.35 

0.05 

0.05 

36 

41 

34 

0.05 

0.35 

0.05 

0.65 

59 

60 

57 

0.20 

0.35 

0.65 

0.05 

22 

23 

22 

0.20 

0.35 

0.50 

0.05 

30 

32 

30 

0.20 

0.35 

0.50 

0.35 

40 

41 

40 

0.20 

0.35 

0.50 

0.20 

56 

58 

56 

0.35 

0.35 

0.05 

0.20 

23 

25 

24 

0.35 

0.35 

0.05 

0.65 

33 

36 

33 

0.35 

0.35 

0.05 

0.50 

37 

41 

37 

0.35 

0.20 

0.65 

0.20 

46 

48 

46 

0.35 

0.50 

0.50 

0.05 

19 

21 

19 

0.35 

0.50 

0.05 

0.50 

28 

31 

28 

0.35 

0.50 

0.05 

0.80 

35 

39 

36 

0.35 

0.50 

0.65 

0.20 

51 

53 

51 

APPENDIX  C 
CHAPTER  3  GRAPHS 


Selected  graphs  of  the  null  power  functions  associated  with  the  comparison  of  two 
independent  multinomial  distributions  are  provided.  Figures  C.l  through  C.3  corre- 
spond to  Fishers  exact  test  with  sample  sizes  30,  40,  and  50  respectively.  Figures  C.4 
through  C.6  represent  the  null  power  using  exact  unconditional  methods  with  sample 
sizes  30,  40,  and  50  respectively. 
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Figure  C.l  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  30. 


Figure  C.2  Null  power  function  for  Fishers  Exact  Test,  for  the 
of  two  independent  multinomial  distributions  when  n  =  40. 
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Figure  C.3  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  50. 
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Figure  C.4  Null  power  function  for  Pearson's  chi-square  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  30. 


83 


Figure  C.5  Null  power  function  for  Pearson's  chi-square  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  40. 
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Figure  C.6  Null  power  function  for  Pearson's  chi-square  Test,  for  the  comparison 
of  two  independent  multinomial  distributions  when  n  =  50. 


APPENDIX  D 
CHAPTER  4  TABLES 


The  following  are  legends  to  be  used  with  the  tables  in  this  appendix. 
Table  D.l 

n     Sample  size  for  each  trinomial  population. 

Ma  Boundary  of  the  critical  region  associated  with  the  exact  unconditional  test. 
a*    Approximate  (to  within  0.001)  size  of  the  test  based  on  Ta. 

Table  D.2 

ne    Sample  size  required  for  the  exact  unconditional  test  statistic. 

Ma    Boundary  of  the  critical  region  associated  with  the  exact  unconditional  test. 

0     Power  for  given  ne  and  parameter  values  0X  and  0y. 

Table  D.3 

nc    Sample  size  required  for  FET  statistic. 

(5     Power  for  given  nc  and  parameter  values  0X  and  6y. 

Table  D.4 

n  Sample  size  for  each  trinomial  population. 

Nn  Number  of  sample  points  for  given  n. 

Nb  Number  of  critical  region  points  in  common. 

Nc  Number  of  points  unique  to  FET  critical  region. 

Ne  Number  of  points  unique  to  exact  unconditional  critical  region. 

Pe  The  proportion  NN<  . 
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Table  D.5 

ne    Sample  size  required  for  the  exact  unconditional  test  statistic. 
nw    Sample  size  required  for  the  asymptotic  test  statistic  with  cc. 
nwo    Sample  size  required  for  the  asymptotic  test  statistic  without  cc. 
nc    Sample  size  required  for  FET  statistic. 
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Table  D.l:  MQ  and  a*  values  for  the  comparison  of  two  independent  binomial 
distributions  at  two  strata. 


a  = 

0.025 

a  =  0.05 

n 

Ma 

a* 

n 

MQ 

a* 

10 

4.76 

0.02455 

10 

3.66 

0.04817 

11 

4.64 

0.02442 

11 

3.69 

0.04828 

12 

5.13 

0.02418 

12 

3.98 

0.04800 

13 

4.93 

0.02439 

13 

3.86 

0.04874 

14 

4.87 

0.02486 

14 

3.84 

0.04956 

15 

5.25 

0.02392 

15 

4.18 

0.04735 

16 

5.07 

0.02480 

16 

4.01 

0.04816 

17 

4.97 

0.02494 

17 

3.97 

0.04930 

18 

5.41 

0.02371 

18 

4.00 

0.04867 

19 

5.22 

0.02456 

19 

4.20 

0.04901 

20 

5.12 

0.02502 

20 

4.12 

0.04930 

21 

5.07 

0.02499 

21 

4.11 

0.04832 

22 

5.40 

0.02481 

22 

4.12 

0.04982 

23 

5.28 

0.02481 

23 

4.29 

0.04949 

24 

5.21 

0.02499 

24 

4.22 

0.04925 

25 

5.19 

0.02501 

25 

4.19 

0.04942 

26 

5.46 

0.02491 

26 

4.15 

0.04978 

27 

5.35 

0.02486 

27 

4.42 

0.04810 

28 

5.30 

0.02493 

28 

4.33 

0.04936 

29 

5.28 

0.02489 

29 

4.28 

0.04907 

30 

5.56 

0.02462 

30 

4.26 

0.04955 

31 

5.44 

0.02483 

31 

4.22 

0.04996 

32 

5.39 

0.02496 

32 

4.45 

0.04920 

33 

5.37 

0.02472 

33 

4.37 

0.05003 

34 

5.69 

0.02433 

34 

4.33 

0.04988 

35 

5.55 

0.02503 

35 

4.31 

0.04992 

36 

5.48 

0.02492 

36 

4.29 

0.04992 

37 

5.45 

0.02489 

37 

4.52 

0.04914 

38 

5.42 

0.02494 

38 

4.44 

0.04926 

Table  D.l  (con't) 


Ui  — 

rv  —  0  05 

M 

M 

lvla 

a* 

Oct 

O.  1  U 

0  09487 

4  39 

0  05001 

40 

4U 

o.oi) 

D  094Q7 

40 

4  36 

0  04936 

41 

5  54 

U .  tJ4 

D  09493 

41 

4  35 

0  05002 

49 

4Zi 

5  5D 

0  09485 

42 

4  36 

0  04938 

4^ 

40 

^  48 

0.40 

0  0940^ 

43 

40 

4  59 

0  04907 

44 

^  7fi 

0  094Q9 

44 

44 

4  47 

4.4  / 

0  04Qfi1 

45 

5.65 

a  no  1  on 

0.02482 

45 

A    A  A 

4.44 

0.04961 

46 

5.59 

0.02502 

46 

4.41 

0.04959 

47 

5.55 

0.02500 

47 

4.40 

0.04961 

48 

5.53 

0.02497 

48 

4.45 

0.04947 

49 

5.53 

0.02484 

49 

4.56 

0.04968 

50 

5.74 

0.02481 

50 

4.51 

0.04971 
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Table  D.2:  Sample  sizes  necessary  to  achieve  a  power  of  at  least  0.80,  at  level 
a  =  0.025,  for  comparing  two  independent  binomial  distributions  at  two  strata, 
using  exact  unconditional  methods. 


01 

09 

Ai 

A2 

rip 

Ma 

a* 

0 

0.75 

0.75 

0.25 

0.25 

11 

4.64 

0.02442 

0.8178 

0.75 

0.75 

0.30 

0.30 

14 

4.87 

0.02486 

0.8299 

0.75 

0.75 

0.35 

0.35 

17 

4.97 

0.02494 

0.8076 

0.75 

0.75 

0.40 

0.40 

22 

5.40 

0.02481 

0.8133 

0.75 

0.75 

0.40 

0.45 

25 

5.19 

0.02501 

0.8030 

0.75 

0.75 

0.45 

0.45 

30 

5.56 

0.02462 

0.8137 

0.75 

0.75 

0.70 

0.30 

31 

5.44 

0.02492 

0.8102 

0.75 

0.75 

0.50 

0.50 

42 

5.50 

0.02485 

0.8143 

0.75 

0.75 

0.70 

0.35 

38 

5.42 

0.02493 

0.8072 

0.75 

0.75 

0.75 

0.35 

40 

5.59 

0.02502 

0.8091 

0.90 

0.80 

0.50 

0.50 

20 

5.12 

0.02502 

0.8317 

0.90 

0.80 

0.60 

0.50 

25 

5.19 

0.02501 

0.8102 

0.90 

0.80 

0.60 

0.55 

30 

5.56 

0.02462 

0.8042 

0.90 

0.80 

0.60 

0.60 

35 

5.55 

0.02503 

0.8045 

0.90 

0.80 

0.65 

0.60 

42 

5.50 

0.02485 

0.8083 

0.95 

0.65 

0.50 

0.35 

17 

4.97 

0.02494 

0.8197 

0.95 

0.65 

0.55 

0.35 

20 

5.12 

0.02502 

0.8320 

0.95 

0.65 

0.65 

0.35 

25 

5.19 

0.02501 

0.8140 

0.95 

0.65 

0.65 

0.40 

30 

5.28 

0.02489 

0.8036 

0.95 

0.65 

0.65 

0.45 

36 

5.48 

0.02492 

0.8169 

0.95 

0.65 

0.65 

0.50 

42 

5.50 

0.02485 

0.8102 
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Table  D.3:  Sample  sizes  necessary  to  achieve  a  power  of  at  least  0.80,  at  level 
a  =  0.025,  for  comparing  two  independent  binomial  distributions  at  two  strata, 
using  conditional  methods. 


u  1 

00 

Ai 

A2 

is 

0 

0.75 

0.75 

0.25 

0.25 

12 

0.8095 

0.75 

0.75 

0.30 

0.30 

15 

0.8323 

0.75 

0.75 

0.35 

0.35 

18 

0.8010 

0.75 

0.75 

0.40 

0.40 

24 

0.8169 

0  75 

0  75 

0.40 

0.45 

27 

0.8106 

0.75 

0.75 

0.45 

0.45 

31 

0.8045 

0.75 

0.75 

0.70 

0.30 

27 

0.8002 

0.75 

0.75 

0.50 

0.50 

44 

0.8065 

0.75 

0.75 

0.70 

0.35 

34 

0.8025 

0.75 

0.75 

0.75 

0.35 

35 

0.8095 

0.90 

0.80 

0.50 

0.50 

20 

0.8116 

0.90 

0.80 

0.60 

0.50 

27 

0.8165 

0.90 

0.80 

0.60 

0.55 

30 

0.8052 

0.90 

0.80 

0.60 

0.60 

34 

0.8063 

0.90 

0.80 

0.65 

0.60 

42 

0.8085 

0.95 

0.65 

0.50 

0.35 

16 

0.8037 

0.95 

0.65 

0.55 

0.35 

19 

0.8293 

0.95 

0.65 

0.65 

0.35 

25 

0.8131 

0.95 

0.65 

0.65 

0.40 

27 

0.8020 

0.95 

0.65 

0.65 

0.45 

30 

0.8007 

0.95 

0.65 

0.65 

0.50 

33 

0.8066 

Table  D.4:  A  comparison  of  the  FET  and  unconditional  critical  regions. 


n 

A  J 

Nn 

Nb 

AT 

M 

t  e 

10 

14641 

6396 

COO 

592 

ooO 

A  Q  ft 

4o.b 

11 

20736 

9816 

T  1  O 

712 

O  1  £ 

816 

ro  A 

53.4 

12 

28561 

-i    A  f~\  o  o 

14088 

1  O  f*  A 

1264 

70/; 

736 

OCt  o 

3b. 8 

13 

38416 

20280 

1  O  T/° 

1376 

1  1  O  A 

1124 

A  C  A 

45.0 

14 

50625 

OO  O  1  o 

28212 

1  7£JA 

1760 

1  /I  OO 

1432 

a  k  n 
40. U 

15 

65536 

O  TO  O  O 

37388 

25b8 

1  O  1  O 

lol2 

OO  Q 
OO.O 

16 

83521 

a  no  1  c 

49816 

Orion 
30o0 

iyzu 

QQ  /I 

oo.4 

i  T 

17 

104976 

O  A  £?0  A 

64624 

oono 
3392 

24U0 

/1 1  /I 
41.4 

18 

130321 

01  coo 

81528 

cc  o  1  o 

5312 

iy40 

Oft  Q 

i  o 

19 

1  O  C\C\C\C\ 

160000 

1  000£?A 

103260 

5560 

2o4U 

OO  Q 
OO.O 

20 

1  O  A  A  O  *t 

194481 

1  OOOO  A 

128924 

6376 

OO  1  o 

3812 

3  ( .4 

O  1 

21 

oo  /for"/? 

234256 

1  c:  O  A  TO 

158472 

7H70 

7072 

4byb 

on  o 

oy.y 

OO 

22 

O  TO  O  /I  1 

279841 

1  O  1  O  1  £} 

191216 

9936 

a  r\c\a 
40yb 

OO  o 

2y.2 

o  o 

23 

O  O  1  TT/? 

331776 

OO  1  f\CO 

231068 

1  A7CO 

10768 

a  nno 

4yy2 

01  1 

ol.  ( 

o 

24 

OOO/^O  c 

390625 

OTZ?  OO  O 

276868 

1  OOOO 

12320 

a  roc 

b53b 

O  1  7 

34.  / 

O  v* 

25 

456976 

OOTO  /I  O 

327943 

13536 

7808 

oft  ft 
3b.b 

26 

531441 

384568 

i  tto ty 

17736 

TOOO 

7000 

OO  O 

28.3 

27 

614656 

451432 

1  o  o  r\  r\ 

18200 

O  /i  OO 

8420 

31. b 

28 

707281 

r"  o/° *i  o 

526612 

1  OOTi° 

19976 

OTOO 

9720 

OO  T 

32.7 

oo 

29 

O  "1  oooo 

810000 

Z1  O  OO  o  o 

608932 

21 352 

1  OOOO 

10988 

O  A  O 

34.0 

30 

ooo  r*  o  i 

923521 

O OTOOO 

697232 

OO  O  A  O 

28640 

OO  A  O 

9040 

O  /I  o 

24.0 

31 

1048576 

OOOO  ,4  o 

800248 

OO  O  OO 

29320 

1  OOOO 

12032 

OO  1 

29.1 

32 

i 1 oc OO 1 

1185921 

C\A  A  OOO 

914932 

O  1  OOO 

31980 

1  A  1  1  O 

14112 

o  a  a 

30. b 

33 

1336336 

1038300 

34448 

15608 

31.2 

34 

1500625 

1668300 

45640 

11508 

20.1 

35 

1679616 

1320408 

45240 

15352 

25.3 

36 

1874161 

1486376 

48088 

18380 

27.5 

37 

2085136 

1663340 

51488 

21176 

29.1 

38 

2313441 

2313441 

55304 

24396 

30.6 

Table  D.4  (con't) 


n 

Nn 

Nb 

Nc 

Ne 

Pe 

39 

2560000 

2560000 

69604 

18816 

21.3 

40 

2825761 

2825761 

70720 

23120 

24.6 

41 

3111696 

3111696 

73560 

26352 

26.4 

42 

3418801 

3418801 

78432 

29632 

27.4 

43 

3748096 

3748096 

83080 

33432 

28.7 

44 

4100625 

4100625 

105520 

25944 

19.7 

45 

4477456 

4477456 

105296 

32940 

23.8 

46 

4879681 

4879681 

109208 

37480 

25.6 

47 

5308416 

5308416 

113352 

41452 

26.8 

48 

5764801 

5764801 

120896 

45120 

27.2 

49 

6250000 

6250000 

127264 

48668 

27.7 

50 

6765201 

6765201 

151752 

40644 

21.1 
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Table  D.5:  Comparisons  between  exact  unconditional,  asymptotic,  and  conditional 
sample  sizes. 


01 

02 

Ai 

A2 

ne 

nc 

nw 

0.75 

0.75 

0.25 

0.25 

11 

12 

10 

12 

0.75 

0.75 

0.30 

0.30 

14 

15 

12 

14 

0.75 

0.75 

0.35 

0.35 

17 

19 

15 

18 

0.75 

0.75 

0.40 

0.40 

22 

24 

19 

22 

0.75 

0.75 

0.40 

0.45 

25 

27 

22 

25 

0.75 

0.75 
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APPENDIX  E 
CHAPTER  4  GRAPHS 


Selected  graphs  of  the  null  power  functions  associated  with  the  comparison  of  two 
independent  binomial  distributions  at  two  different  strata  are  provided.  Figures  E.l 
through  E.5  correspond  to  Fishers  exact  test  with  sample  sizes  10,  20,  30,  40,  and  50 
respectively.  Figures  E.6  through  E.10  represent  the  null  power  using  exact  uncondi- 
tional methods  with  sample  sizes  10,  20,  30,  40,  and  50  respectively. 
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Figure  E.l  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  binomial  distributions  at  two  strata  when  n  =  10. 
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Figure  E.2  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  binomial  distributions  at  two  strata  when  n  =  20. 


Figure  E.3  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  binomial  distributions  at  two  strata  when  n  =  30. 


Figure  E.4  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  binomial  distributions  at  two  strata  when  n  =  40. 


Figure  E.5  Null  power  function  for  Fishers  Exact  Test,  for  the  comparison 
of  two  independent  binomial  distributions  at  two  strata  when  n  =  50. 
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Figure  E.6  Null  power  function  for  MH  with  continuity  correction,  for  the 
comparison  of  two  independent  binomial  distributions  at  two  strata  when  n 


Figure  E.7  Null  power  function  for  MH  with  continuity  correction,  for  the 
comparison  of  two  independent  binomial  distributions  at  two  strata  when  n 


Figure  E.8  Null  power  function  for  MH  with  continuity  correction,  for  the 
comparison  of  two  independent  binomial  distributions  at  two  strata  when  n 
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Figure  E.9  Null  power  function  for  MH  with  continuity  correction,  for  the 
comparison  of  two  independent  binomial  distributions  at  two  strata  when  n  =  40. 


Figure  E.10  Null  power  function  for  MH  with  continuity  correction,  for  the 
comparison  of  two  independent  binomial  distributions  at  two  strata  when  n  = 


APPENDIX  F 
FORTRAN  PROGRAMS 


The  following  Fortran  programs  are  provided  here  : 

1.  bnd2by3  :  computes  the  exact  unconditional  critical  points  for  the  2  x  3 
contingency  table. 

2.  fet2by3crit  :  determines  the  critical  regions  for  FET  for  the  2  x  3 
contingency  table. 

3.  fet2by3pwr  :  computes  power  for  FET  for  the  2  x  3  case. 

4.  bnd2by2by2  :  computes  the  exact  unconditional  critical  points  for  the  2  x  2  x  2 
contingency  table. 

5.  fet2by2by2crit  :  determines  the  critical  regions  for  FET  for  the  2  x  2  x  2 
contingency  table. 

6.  fet2by2by2pwr  :  computes  power  for  FET  for  the  2  x  2  x  2  case. 
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C 


program  bnd2by3 

integer  mm,mmm,critcnt,i,c(100),m,cnt,totcnt,j,k,ii,jj 
double  precision  xl,x2,x3,yl,y2,y3,n,bnd,nx,ny,cl,c2 
double  precision  prbl,prb2,pbndl,pbnd2,pbnd3,n2,c3,c4 
double  precision  c5,c6,c7,bndprb,xlow,xhigh,ylow,yhigh 
double  precision  delta(1000),xxl(1000000),xx2(1000000) 
double  precision  yyl(1000000),yy2(1000000),bound(1000000) 
double  precision  lowl(10,100000),low2(10,100000),hl 
double  precisionpl,p2,p3,h2,a,b,c62,c8,bndprb2 
double  precision  bound2(100000),mxbnd,mxfnct 

read*,  n,bnd 
cnt  =  1 
k  =  0 

n2  =  2.0D0*n 

delta(l)  =  1.0D0/(2.0D0*n2) 
critcnt  =  0 

C 

C    Determine  the  critical  region 

do  xl  =  0.0D0,n 
do  x2  =  0.0D0,n-xl 
do  yl  =  0.0D0,n 
do  y2  =  0.0D0,n-yl 
if  (xl  .EQ.  yl  .AND.  x2  .GE.  y2)  then 
goto  10 
endif 

if  (xl  .LT.  yl)  then 
goto  10 
endif 

x3  =  n  -  xl  -  x2 
y3  =  n  -  yl  -  y2 
tl  =  xl  +  yl 
t2  =  x2  +  y2 
t3  =  x3  4-  y3 
tlsq  =  (xl  -  yl)**2 
t2sq  =  (x2  -  y2)**2 
t3sq  =  (x3  -  y3)**2 

if  (tl  .EQ.  0.0  .OR.  t2  .EQ.  0.0  .OR.  t3  .EQ.  0.0)  then 
ctot  =  tlsq/(tl+0.001)  +  t2sq/(t2+0.001)  +  t3sq/(t3+0.001) 

plfiP 

ctot  =  tlsq/tl  +  t2sq/t2  +  t3sq/t3 
endif 

if  (ctot  .GE.  bnd)  then 
critcnt  =  critcnt  +  1 
xxl  (critcnt)  =  xl 
xx2  (critcnt)  =  x2 
yyl  (critcnt)  =  yl 
yy2(critcnt)  =  y2 


endif 
10  end  do 
end  do 
end  do 
end  do 

Q 

C  First  attempt  to  bound  function 
C 

mmm  =  0 
totcnt  =  0 
mm  =  0 
m  =  1 
c(l)  =  0 

hi  =  0.50D0  -  delta(m) 
do  xlow  =  0.0D0,hl,delta(m) 
xhigh  =  xlow  +  delta(m) 
pi  =  xlow  +  delta(m)/2.0D0 
do  ylow  =  0.0D0,xlow,delta(m) 
yhigh  =  ylow  +  delta(m) 
k  =  k  +i 

p2  =  ylow  +  delta(m)/2.0D0 
p3  =  1.0D0  -  pi  -  p2 
totcnt  =  totcnt  +  1 
bndprb  =  0.0D0 
bndprb2  =  0.0D0 
do  i  =  l,critcnt 
tl  =  xxl(i)+yyl(i) 
t2  =  xx2(i)+yy2(i) 

t3  =  n2  -  tl  -  t2 
nx  =  n  -  xxl(i) 

ny  =  n  -  yyl(i) 

prbl  =  tl/n2 

prb2  =  t2/n2 

if  (prbl  .LE.  xlow)  then 

pbndl  -  xlow 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 

endif 
goto  20 

else 

pbndl  =  xhigh 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 


endif 

endif  ,       ,  ,„ 

20  pbnd3  =  1.0D0  -  pbndl  -  pbnd2 

call  combo(n,xxl(i),cl) 

call  combo(nx,xx2(i),c2) 

call  combo(n,yyl(i),c3) 

call  combo(ny,yy2(i),c4) 

c5  _  cl*c2*c3*c4 

c6  =  (Pbndl**tl)*(pbnd2**t2)*(pbnd3**t3) 
c62  =  (Pl**tl)*(p2**t2)*(p3**t3) 

c7  =  c5*c6 

c8  =  c5*c62 

bndprb  =  bndprb  +  c7 

bndprb2  =  bndprb2  +  c8 

end  do 

bndprb  =  2.0D0*bndprb 
bound2(k)  =  2.0D0*bndprb2 

Q 

C    Determine  which  regions  are  okay  and  which  need 

C    to  be  partitioned  further 

C 

if  (bndprb  .LE.  0.0250D0)  then 
mmm  =  mmm  +  1 
bound(mmm)  =  bndprb 
else 

c(cnt)  =  c(cnt)  +  1 

lowl(m,c(cnt))  =  xlow 

low2(m,c(cnt))  =  ylow 

endif 
end  do 
end  do 

if  (c(cnt)  .GT.  0)  then 
m  =  m  +  1 

delta(m)  =  delta(m-l)/2.0D0 

goto  30 

endif 

C 

C    Determining  bounds  for  the  above  regions  (  that  needed 
C    further  partitioning  )  iteratively,  until  the  desired 
C    accuracy  is  attained. 
C 

30  c(cnt+l)  =  0 
do  mm  =  l,c(cnt) 
xxlow(l)  =  lowl(m-l,mm) 
xxlow(2)  =  lowl(m-l,mm)  +  delta(m-l) 
yylow(l)  =  low2(m-l,mm) 
yylow(2)  =  low2(m-l,mm)  +  delta(m-l) 
do  ii  =  1,2 


do  jj  =  1,2 

xlow  =  xxlow(ii) 

xhigh  =  xlow  +  delta(m) 

pi  =  xlow  +  delta(m)/2.0D0 

ylow  =  yylowQj) 

yhigh  =  ylow  +  delta(m) 

p2  =  ylow  +  delta(m)/2.0D0 

p3  =  1.0D0  -  pi  -  p2 

k  =  k  +  1 

bndprb  =  0.0D0 

bndprb2  =  0.0D0 

do  i  =  l,critcnt 

tl  =  xxl(i)+yyl(i) 

t2  =  xx2(i)+yy2(i) 

t3  =  n2  -  tl  - 12 
nx  =  n  -  xxl(i) 

ny  =  n  -  yyl(i) 

prbl  =  tl/n2 

prb2  =  t2/n2 

if  (prbl  .LE.  xlow)  then 

pbndl  =  xlow 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 

endif 
goto  25 

else 

pbndl  =  xhigh 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 

endif 

endif 

25  pbnd3  =  1.0D0  -  pbndl  -  pbnd2 
call  combo(n,xxl(i),cl) 
call  combo(nx,xx2(i),c2) 
call  combo(n,yyl(i),c3) 
call  combo(ny,yy2(i),c4) 
c5  _  cl*c2*c3*c4 

c6  =  (Pbndl**tl)*(pbnd2**t2)*(pbnd3* 

c62  =  (pl**tl)*(p2**t2)*(p3**t3) 

c7  —  c5*c6 

c8  =  c5*c62 

bndprb  =  bndprb  +  c7 

bndprb2  =  bndprb2  +  c8 
end  do 


bndprb  =  2.0D0*bndprb 
bound2(k)  =  2.0D0*bndprb2 
if  (bndprb  .LE.  0.0250D0)  then 
mmm  =  mmm  +  1 
bound(mmm)  =  bndprb 
else 

c(cnt+l)  =  c(cnt+l)  +  1 

lowl(m,c(cnt+l))  =  xlow 

low2(m,c(cnt+l))  =  ylow 

endif 
end  do 
end  do 
end  do 

cnt  =  cnt  +  1 

if  (c(cnt)  .GT.  0)  then 
m  =  m  +  1 

delta(m)  =  delta(m-l)/2.0D0 

goto  30 

endif 

C 

C    Determining  accuracy  of  the  bound. 
C 

40  a  =  bound(l) 
do  j  =  2, mmm 
if  (bound(j)  .GT.  a)  then 
b  =  bound  (j) 

endif 

a  =  b 
end  do 

mxbnd  =  a 
a  =  bound2(l) 

do  j  =  2,k 

if  (bound2(j)  .GT.  a)  then 
b  =  bound2(j) 

endif 

a  =  b 
end  do 

mxfnct  =  a 

print*, mxbnd, mxbnd-mxfnct 
end 

subroutine  combo  (nn,rr, choose) 
double  precision  nn,rr,cc,diff,choose,dd, 
diff  =  nn  -  rr 

if  (rr  .EQ.  0.0D0  .OR.  rr  .EQ.  nn)  then 
choose  =  1.0D0 

6  Is  6 

if  (diff  .GE.  rr)  then 
dd  =  1.0D0 


cc  =  l.ODO 

do  i  =  nn,diff+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,rr 

dd  =  dd*i 
end  do 

choose  =  cc/dd 

else 

dd  =  l.ODO 
cc  =  l.ODO 

do  i  =  nn,rr+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,diff 

dd  =  dd*i 
end  do 

choose  =  cc/dd 

endif 
endif 
return 
end 
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program  fet2by3crit 
character  *9  fill,fil2 
character  *3  fisher, numbr 
ch.3.r3.ctGr  ^4  extn 

integer  f,j,jj,i,o,ff,nmber(900000),count,cl,c2,c3,cntl(900000) 

integer  cnt2(900000),cnt3(10000),indx(900000),k,l,ntemp 

double  precision  suml,sum2,sum3,xl,x2,x3,yl,y2,y3,n,s,flag,tl,t2 

double  precision  rx,ry,a,b,c,e,n2,nx,ny,nt,s2(900000),s3(900000) 

double  precision  diffl,diff2,diff3,d,perm(900000),ttl,tt2 

double  precision  rl(900000),r2(900000),r3(900000),sl(900000),t3 

double  precision  xxl,xx2,yyl,yy2,g,h,q,tt,axl,ax2,ax3,ayl,ay2,ay3 

double  precision  numl,num2,num3,num4,num,denl 

double  precision  den2,den,probl  (900000), alpha,prbtype(900000) 

alpha  =  0.025D0 

do  k  =  1,6 

do  1  =  0,9 

fisher  =  'fet' 

numbr  =  'nmb' 

extn  =  '.dat' 

fill  =  fisher  //  char(48+k)  //  char (48+1)  //  extn 
fil2  =  numbr  //  char(48+k)  //  char (48+1)  //extn 
print*, fill 

ntemp  =  (10*k)  +  1 

n  =  dble(ntemp) 

n2  =  2.0D0*n 

tot  =  0.0D0 

count  =  0 

cl  =  0 

c2  =  0 

c3  =  0 

tt  =  0.0D0 

q  =  int(n/2.0D0) 

ff  =  0 
j  =  0 
o  =  0 
f  =  0 

C 

C    Determine  Un,  and  the  representations  for  each  element  of  Un 
C    as  was  described  in  Chapter  3. 


C 


do  xxl  =  0.0D0,q 
do  xx2  =  0.0D0,n-xxl 
do  yvl  =  0.0D0,q 
do  yy2  =  0.0D0,n-yyl 
flag  =  0.0D0 
jj  =  0 
xl  =  xxl 
x2  =  xx2 


yi  =  yyi 
y2  =  yy2 

x3  =  n  -  xl  -  x2 
y3  =  n  -  yl  -  y2 

rx  =  max(xl,x2,x3) 

ry  =  max(yl,y2,y3) 

if  (rx  .It.  ry)  then 

a  =  xl 
b  =  x2 
c  =  x3 
e  =  yl 

g  =  y2 
h  =  y3 


x3  =  h 
yl  =  a 

y2  =  b 

y3  =  c 

endif 

if  (rx  .eq.  ry)  then 

axl  =  xl 
ax2  =  x2 
ax3  =  x3 
ayl  =  yl 
ay2  =  y2 
ay3  =  y3 

if  (axl  .It.  ax3)  then 

a  =  axl 

b  =  ax3 

axl  =  b 

ax3  -  a 

endif 

if  (axl  .It.  ax2)  then 

a  =  axl 
b  =  ax2 
axl  =  b 
ax  2  =  a 
endif 

if  (ax2  .It.  ax3)  then 

a  =  ax2 

b  =  ax3 

ax2  =  b 

ax3  =  a 

endif 

if  (ayl  .It.  ay3)  then 
a  =  ayl 
b  =  ay3 


ayl  =  b 
ay3  -  a 
endif 

if  (ayl  .It.  ay2)  then 

a  =  ayl 

b  =  ay2 

ayl  =  b 

ay2  =  a 

endif 

if  (ay2  .It.  ay3)  then 

a  =  ay2 

b  =  ay3 

ay2  =  b 

ay3  =  a 

endif 

if  (ay2  .It.  ax2)  then 
a  =  xl 
b  =  x2 
c  =  x3 
e  =  yl 

g  =  y2 

h  =  y3 
xl  =  e 
x2  =  g 

x3  =  h 
yl  =  a 

v2  =  b 

y3  =  c 

endif 

if  (ax2  .eq.  ay2)  then 

flag  =  l.ODO 

endif 
endif 

20  if  (xl  .It.  x3)  then 
a  =  xl 
b  =  x3 
c  =  yl 
e  =  y3 

xl  =  b 
x3  —  a 
yl  =  e 

y3  =  c 

endif 

if  (xl  .It.  x2)  then 

a  =  xl 
b  =  x2 
c  =  yl 

e  =  y2 


xl  =  b 
x2  =  a 
yl  =  e 

y2  =  c 

endif 

if  (x2  .It.  x3)  then 
a  =  x2 
b  =  x3 
c  =  y2 
e  =  y3 

x2  =  b 
x3  =  a 
y2  =  e 

y3  =  c 

endif 

if  (xl  .eq.  x2  .and.  yl  .It.  y2)  then 

a  =  yl 

b  =  y2 

yl  =  b 

y2  =  a 

endif 

if  (x2  .eq.  x3  .and.  y2  .It.  y3)  then 

a  =  y2 

b  =  y3 

y2  =  b 

y3  =  a 

endif 

if  (flag  .eq.  l.ODO)  then 
if  (xl  .ne.  yl)  then 
tolx  =  abs(xl  -  yl) 
if  (y2  .eq.  ry)  then 
toly  =  abs(y2  -x2) 
endif 

if  (y3  .eq.  ry)  then 

toly  =  abs(y3  -  x3) 

endif 
endif 

if  (tolx  .It.  toly)  then 

a  =  xl 
b  =  x2 
c  =  x3 
e  =  yl 

g  =  y2 
h  =  y3 

xl  =  e 
x2  =  g 

x3  =  h 
yl  =  a 
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y2  =  b 
y3  =  c 
goto  20 
endif 
endif 

if  (f  .eq.  0)  then 

j  =  j  +  1 
rl(j)  =  xl 
r2(j)  =  x2 
r3(j)  =  x3 
sl(j)  =  yi 
s20)  =  y2 
s3(j)  =  y3 
goto  10 
endif 
do  i  =  l,j 

if  (xl  .eq.  rl(i)  .and.  x2  .eq.  r2(i))  then 
if  (yl  .eq.  si (i)  .and.  y2  .eq.  s2(i))  then 

jj  =  jj  +  1 
endif 
endif 
end  do 

if  ( jj  .eq.  0)  then 

j  =  j  +  1 
rl(j)  =  xl 
r2(j)  =  x2 
r3(j)  =  x3 

sl(j)  =  yi 
s2(j)  =  y2 
s3(j)  =  y3 
end  if 
10  f  =  f  +  1 
end  do 
end  do 
end  do 
end  do 
do  i  =  l,j 
flag  =  0.0D0 
suml  =  rl(i)  +  sl(i) 
sum2  =  r2(i)  +  s2(i) 
sum3  =  r3(i)  +  s3(i) 
diffl  =  abs(rl(i)  -  sl(i)) 
diff2  =  abs(r2(i)  -  s2(i)j 
diff3  =  abs(r3(i)  -  s3(i)) 
if  (rl(i)  .eq.  s2(i)  .and.  r2(i)  .eq.  sl(i))  then 
flag  =  1.0D0 
endif 

if  (rl (i)  .eq.  s3(i)  .and.  r3(i)  .eq.  sl(i))  then 


flag  =  l.ODO 

endif  , ,  , . 

if  (r2(i)  .eq.  s3(i)  .and.  r3(i)  .eq.  s2(i))  then 

flag  =  l.ODO 
endif 

if  (r2(i)  .eq.  r3(i)  .and.  s2(i)  .eq.  s3(i))  then 

flag  =  l.ODO 

endif 

if  (rl(i)  .eq.  r3(i)  .and.  sl(i)  .eq.  s3(i))  then 

flag  =  l.ODO 

endif 

if  (rl(i)  .eq.  r2(i)  .and.  sl(i)  .eq.  s2(i))  then 

flag  =  l.ODO 

endif 

if  (suml  .eq.  sum2  .and.  sum2  .eq.  sum3)  then 

s  =  3.0DO 
endif 

if  (suml  .eq.  sum2  .and.  sum2  .ne.  sum3)  then 

s  =  2.0D0 
endif 

if  (suml  .eq.  sum3  .and.  sum2  .ne.  sum3)  then 

s  =  2.0D0 
endif 

if  (sum2  .eq.  sum3  .and.  suml  .ne.  sum3)  then 

s  =  2.0D0 
endif 

if  (suml  .ne.  sum2  .and.  sum2  .ne.  sum3)  then 

if  (suml  .ne.  sum3)  then 

s  =  l.ODO 

endif 

endif 

if  (diffl  .eq.  diff2  .and.  diff2  .eq.  diff3)  then 

d  =  3.0D0 
endif 

if  (diffl  .eq.  diff2  .and.  diff2  .ne.  diff3)  then 

d  =  2.0D0 
endif 

if  (diffl  .eq.  diff3  .and.  diff2  .ne.  diff3)  then 

d  =  2.0D0 
endif 

if  (diff2  .eq.  diff3  .and.  diffl  .ne.  diff3)  then 

d  =  2.0D0 
endif 

if  (diffl  .ne.  diff2  .and.  diff2  .ne.  diff3)  then 

if  (diffl  .ne.  diff3)  then 

d  =  l.ODO 

endif 

endif 


if  (s  .eq.  l.ODO  .and.  d  .eq.  3.0D0)  then 

perm(i)  =  6.0D0 

endif 

if  (s  .eq.  l.ODO  .and.  d  .eq.  2.0D0)  then 

perm(i)  =  12.0D0 

endif 

if  (s  .eq.  l.ODO  .and.  d  .eq.  l.ODO)  then 

perm(i)  =  12.0D0 

endif 

if  (s  .eq.  3.0D0  .and.  d  .eq.  3.0D0)  then 

perm(i)  =  l.ODO 

endif 

if  (s  .eq.  3.0D0  .and.  d  .eq.  2.0D0)  then 

perm(i)  =  6.0D0 

endif 

if  (s  .eq.  3.0D0  .and.  d  .eq.  l.ODO)  then 

perm(i)  =  12.0D0 

endif 

if  (s  .eq.  2.0D0  .and.  d  .eq.  2.0D0  .and.  flag  .eq.  O.ODO)  then 

perm(i)  =  12.0D0 

endif 

if  (s  .eq.  2.0D0  .and.  d  .eq.  2.0D0  .and.  flag  .eq.  l.ODO)  then 

perm(i)  =  6.0D0 

endif 

if  (s  .eq.  2.0D0  .and.  d  .eq.  l.ODO)  then 

perm(i)  =  12.0D0 

endif 

if  (s  .eq.  2.0DO  .and.  d  .eq.  3.0D0)  then 

perm(i)  =  3.0D0 

endif 
end  do 

do  i  =  1  j 

tl  =  rl(i)  +  sl(i) 

t2  =  r2(i)  +  s2(i) 

t3  =  (2.0D0*n)  -  (tl  +  t2) 

if  (tl  .LT.  t2)  then 

a  =  rl(i) 

b  =  sl(i) 

rl(i)  =  r2(i) 

sl(i)  =  s2(i) 

r2(i)  =  a 

s2(i)  =  b 

tl  =  rl(i)  +  sl(i) 

t2  =  r2(i)  +  s2(i) 

endif 

if  (tl  .LT.  t3)  then 
a  =  rl(i) 
b  =  sl(i) 
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rl(i)  =  r3(i) 

sl(i)  =  s3(i) 

r3(i)  =  a 

s3(i)  =  b 

tl  =  rl(i)  +  sl(i) 

t3  =  r3(i)  +  s3(i) 

endif 

if  (t2  .LT.  t3)  then 

a  =  r2(i) 

b  =  s2(i) 

r2(i)  =  r3(i) 

s2(i)  =  s3(i) 

r3(i)  =  a 

s3(i)  =  b 

t2  =  rl(i)  +  sl(i) 

t3  =  r3(i)  +  s3(i) 

endif 

tl  =  rl(i)  +  sl(i) 
t2  =  r2(i)  +  s2(i) 
t3  =  r3(i)  +  s3(i) 

if  (tl  .EQ.  t2  .AND.  t2  .EQ.  t3)  then 

nmber(i)  =  3 

endif 

if  (tl  .EQ.  t2  .AND.  t2  .NE.  t3)  then 

nmber(i)  =  2 

endif 

if  (t2  .EQ.  t3  .AND.  t2  .NE.  tl)  then 

nmber(i)  =  2 

endif 

if  (tl  .NE.  t2  .AND.  tl  .NE.  t3  .AND.  t2  .NE.  t3)  then 

nmber(i)  =  1 

endif 
end  do 

C 

C    Determine  critical  region 
C 

do  i  =  l,j 
tl  =  rl(i)  +  sl(i) 
t2  =  r2(i)  +  s2(i) 
nx  =  n  -  rl(i) 
ny  =  n  -  sl(i) 
nt  =  n2  -  tl 

call  combo(n,rl(i),numl) 
call  combo(nx,r2(i),num2) 
call  combo(n,sl(i),num3) 
call  combo(ny,s2(i),num4) 
call  combo(n2,tl,denl) 
call  combo(nt,t2,den2) 


num  =  numl*num2*num3*num4 
den  =  denl*den2 
probl(i)  =  num/den 

end  do 

do  i  =  l,j 

tl  =  rl(i)  +  sl(i) 
t2  =  r2(i)  +  s2(i) 
t3  =  r3(i)  +  s3(i) 

if  (nmber(i)  .EQ.  3  .AND.  perm(i)  .EQ.  1.0D0)  then 
probl(i)  =  probl(i) 

endif  /%  .  . 

if  (nmber(i)  .EQ.  3  .AND.  perm(i)  .EQ.  6.0D0)  then 

probl(i)  =  6.0D0*probl(i) 

endif 

if  (nmber(i)  .EQ.  3  .AND.  perm(i)  .EQ.  12.0D0)  then 

probl(i)  =  12.0D0*probl(i) 

endif 

if  (nmber(i)  .EQ.  2  .AND.  perm(i)  .EQ.  6.0D0)  then 

probl(i)  =  2.0D0*probl(i) 

endif 

if  (nmber(i)  .EQ.  2  .AND.  perm(i)  .EQ.  12.0D0)  then 

probl(i)  =  4.0D0*probl(i) 

endif 

if  (nmber(i)  .EQ.  1  .AND.  perm(i)  .EQ.  12.0D0)  then 

probl(i)  =  2.0D0*probl(i) 

endif 
end  do 

do  i  =  l,j 

tot  =  tot  +  perm(i) 

end  do 

print*, tot 
do  i  =  l,j 

if  (probl(i)  .LE.  alpha)  then 
count  =  count  +  1 
rl(count)  =  rl(i) 
r2  (count)  =  r2(i) 
r3(count)  =  r3(i) 
si  (count)  =  sl(i) 
s2(count)  =  s2(i) 
s3(count)  =  s3(i) 
probl  (count)  =  probl(i) 
nmber(count)  =  nmber(i) 
perm(count)  =  perm(i) 

endif 
end  do 

print*,count 
do  i  =  1, count 
tl  =  rl(i)  +  sl(i) 
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t2  =  r2(i)  +  s2(i) 

t3  =  r3(i)  +  s3(i) 

if  (nmber(i)  .EQ.  1)  then 

cl  =  cl  +  1 

cntl(cl)  =  i 

endif 

if  (nmber(i)  .EQ.  2)  then 
c2  =  c2  +  1 
cnt2(c2)  =  i 
endif 

if  (nmber(i)  .EQ.  3)  then 
c3  =  c3  +  1 
cnt3(c3)  =  i 

endif 
end  do 

print*, cl,c2,c3 
count  =  0 
do  i  =  l,cl 
prb  =  0.0D0 

tl  =  rl(cntl(i))  +  sl(cntl(i)) 
t2  =  r2(cntl(i))  +  s2(cntl(i)) 

if  (tl  .EQ.  9.0  .AND.  t2  .EQ.  6.0  .AND.  n  .EQ.  10.0)  then 

print*  tl,t2,probl(cntl(i)) 

endif 

do  ii  =  l,cl 

ttl  =  rl(cntl(ii))  +  sl(cntl(ii)) 

tt2  =  r2(cntl(ii))  +  s2(cntl(ii)) 

if  (tl  .EQ.  ttl  .AND.  t2  .EQ.  tt2)  then 

if  (probl(cntl(ii))  .LE.  probl(cntl(i)))  then 

prb  =  prb  +  probl(cntl(ii)) 

endif 
endif 
end  do 

if  (prb  .LE.  alpha)  then 
count  =  count  +  1 
indx(count)  =  cntl(i) 

endif 
end  do 

do  i  =  l,c2 
prb  =  0.0D0 

tl  =  rl(cnt2(i))  +  sl(cnt2(i)) 
t2  =  r2(cnt2(i))  +  s2(cnt2(i)) 
do  ii  =  l,c2 

ttl  =  rl(cnt2(ii))  +  sl(cnt2(ii)) 

tt2  =  r2(cnt2(ii))  4-  s2(cnt2(ii)) 

if  (tl  .EQ.  ttl  .AND.  t2  .EQ.  tt2)  then 

if  (probl(cnt2(ii))  .LE.  probl(cnt2(i)))  then 

prb  —  prb  +  probl(cnt2(ii)) 


endif 
endif 
end  do 

if  (prb  .LE.  alpha)  then 
count  =  count  +  1 
indx(count)  =  cnt2(i) 

endif 
end  do 

do  i  =  l,c3 
prb  =  0.0D0 

tl  =  rl(cnt3(i))  +  sl(cnt3(i)) 
t2  =  r2(cnt3(i))  +  s2(cnt3(i)) 
do  ii  =  l,c3 

ttl  =  rl(cnt3(ii))  +  sl(cnt3(ii)) 

tt2  =  r2(cnt3(iij)  +  s2(cnt3(ii)) 

if  (tl  .EQ.  ttl  .AND.  t2  .EQ.  tt2)  then 

if  (probl(cnt3(ii))  .LE.  probl(cnt3(i)))  then 

prb  =  prb  +  probl(cnt3(ii)) 

endif 
endif 
end  do 

if  (prb  .LE.  alpha)  then 
count  =  count  +  1 
indx(count)  =  cnt3(i) 

endif 
end  do 

do  i  =  1, count 

tl  =  rl(indx(i))  +  sl(indx(i)) 
t2  =  r2(indx(i))  +  s2(indx(i)) 
t3  =  n2  -  tl  - 12 

if  (perm(indx(i))  .EQ.  12.0D0)  then 

prbtype(i)  =  12.0D0 

endif 

if  (perm(indx(i))  .EQ.  6.0D0  .AND.  tl  .EQ.  t2)  then 

if  (rl(indx(i))  .EQ.  r2(indx(i))  .AND.  sl(indx(i))  .EQ.  s2(indx(i)))  then 

prbtype(i)  =  62.0D0 
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prbtype(i)  =  61.0D0 

endif 
endif 

if  (perm(indx(i))  .EQ.  6.0D0  .AND.  t2  .EQ.  t3)  then 

if  (r3(indx(i))  .EQ.  r2(indx(i))  .AND.  s3(indx(i))  .EQ.  s2(indx(i)))  then 

prbtype(i)  =  63.0D0 

else 

prbtype(i)  =  61.0D0 

endif 
endif 

open(l,FILE=fill) 
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write(l,'(5fl2.5)')  rl(indx(i)),r2(indx(i)),sl(indx(i)),s2(indx(i)),prbtype(i) 
end  do 

open(l,FILE=fil2) 

write(l,'(ilO)')  count 

end  do 
end  do 
end 

subroutine  combo  (nn,rr, choose) 
double  precision  nn,rr,cc,diff,choose,dd,i 
diff  =  nn  -  rr 

if  (rr  .EQ.  0.0D0  .OR.  rr  .EQ.  nn)  then 

choose  =  1.0D0 
else 

if  (diff  .GE.  rr)  then 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,diff+1.0DO,-1.0DO 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,rr 

dd  =  dd*i 
end  do 

choose  =  cc/dd 
else 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,rr+l.OD0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,diff 

dd  =  dd*i 
end  do 

choose  =  cc/dd 

endif 
endif 
return 
end 


program  fet2by3pwr 
integer  m,c,count,j 

double  precision  xl,x2,yl,y2,perm,thetaxl,thetax2,d,e 
double  precision  thetax3,thetay3,n,nx,ny,cxl,cx2,cyl 
double  precision  al,a2,a3,a4,a5,a6,x3,y3,cy2,fctl,fct3 
double  precision  cl,c2,c3,c4,c5,c6,prb,powr,thetayl,thetay2 
read*  ,n,thetaxl  ,thetax2,thetay  1  ,thetay2 
powr  =  0.0D0 

m=-°0 

open(  1  ,file= '/usr  /  users /gcrans /shuster/ fortrn/nmb57.dat ' ) 

read(l,'(i50)')  count 

thetax3  =  1.0D0  -  thetaxl  -  thetax2 

thetay3  =  1.0D0  -  thetayl  -  thetay2 

open(l,file=7usr/users/gcrans/shuster/fortrn/fet57.dat') 

do  j  =  1, count 

c  =  c  +  1 

read(l,'(5fl2.5)')  xl,x2,yl,y2,perm 

nx  =  n  -  xl 
ny  =  n  -  yl 
x3  =  nx  -  x2 
y3  =  ny  -  y2 

call  combo(n,xl,cxl) 
call  combo(nx,x2,cx2) 
call  combo(n,yl,cyl) 
call  combo(ny,y2,cy2) 

fctl  =  cxl*cx2*cyl*cy2 
if  (perm  .EQ.  62.0D0)  then 

d  =  fctl*(thetaxl**xl)*(thetax2**x2)*(thetax3**x3) 
e  =  (thetayl**yl)*(thetay2**y2)*(thetay3**y3) 
al  =  d*e 

d  =  fctl*(thetaxl**x3)*(thetax2**x2)*(thetax3**xl) 

e  =  (thetayl**y3)*(thetay2**y2)*(thetay3**yl) 
a2  _  cj*e 

d  =  fctl*(thetaxl**xl)*(thetax2**x3)*(thetax3**x2) 
e  =  (thetayl**yl)*(thetay2**y3)*(thetay3**y2) 
a3  =  d*e 

d  =  fctl*(thetayl**xl)*(thetay2**x2)*(thetay3**x3) 
e  =  (thetaxl**yl)*(thetax2**y2)*(thetax3**y3) 
a4  =  d*e 

d  =  fctl*(thetayl**x3)*(thetay2**x2)*(thetay3**xl) 
e  =  (thetaxl**y3)*(thetax2**y2)*(thetax3**yl) 

d  =  fctl*(thetayl**xl)*(thetay2**x3)*(thetay3**x2) 
e  =  (thetaxl**yl)*(thetax2**y3)*(thetax3**y2) 
a6  =  d*e 

prb  =  al+a2+a3+a4+a5+a6 
endif 


if  (perm  .EQ.  63.0D0)  then 

d  =  fctl*(thetaxl**xl)*(thetax2**x2)*(thetax3**x3) 
e  =  (thetayl**yl)*(thetay2**y2)*(thetay3**y3) 
al  =  d*e 

d  =  fctl*(thetaxl**x3)*(thetax2**x2)*(thetax3**xl) 
e  =  (thetayl**y3)*(thetay2**y2)*(thetay3**yl) 
a2  _  (\*e 

d  =  fctl*(thetaxl**x2)*(thetax2**xl)*(thetax3**x3) 
e  =  (thetayl**y2)*(thetay2**yl)*(thetay3**y3) 
a3  =  d*e 

d  =  fctl*(thetayl**xl)*(thetay2**x2)*(thetay3**x3) 

e  =  (thetaxl**yl)*(thetax2**y2)*(thetax3**y3) 
a4  =  d*e 

d  =  fctl*(thetayl**x3)*(thetay2**x2)*(thetay3**xl) 
e  =  (thetaxl**y3)*(thetax2**y2)*(thetax3**yl) 

d  =  fctl*(thetayl**x2)*(thetay2**xl)*(thetay3**x3) 
e  =  (thetaxl**y2)*(thetax2**yl)*(thetax3**y3) 
a6  =  d*e 

prb  =  al+a2+a3+a4+a5+a6 
endif 

if  (perm  .EQ.  12.0D0  .OR.  perm  .EQ.  61.0D0)  then 
d  =  fctl*(thetaxl**xl)*(thetax2**x2)*(thetax3**x3) 
e  =  (thetayl**yl)*(thetay2**y2)*(thetay3**y3) 
al  =  d*e 

d  =  fctl*(thetaxl**x2)*(thetax2**xl)*(thetax3**x3) 
e  =  (thetayl**y2)*(thetay2**yl)*(thetay3**y3) 
a2  =  d*e 

d  =  fctl*(thetaxl**x2)*(thetax2**x3)*(thetax3**xl) 
e  =  (thetayl**y2)*(thetay2**y3)*(thetay3**yl) 
a3  =  d*e 

d  =  fctl*(thetaxl**x3)*(thetax2**x2)*(thetax3**xl) 
e  =  (thetayl**y3)*(thetay2**y2)*(thetay3**yl) 
a4  =  d*e 

d  =  fctl*(thetaxl**xl)*(thetax2**x3)*(thetax3**x2) 
e  =  (thetayl**yl)*(thetay2**y3)*(thetay3**y2) 

d  =  fctl*(thetaxl**x3)*(thetax2**xl)*(thetax3**x2) 

e  =  (thetayl**y3)*(thetay2**yl)*(thetay3**y2) 
a6  =  d*e 

prb  =  al+a2-(-a3+a4+a5+a6 
endif 

if  (perm  .EQ.  12.0D0)  then 
m  =  m  +  1 

d  =  fctl*(thetayl**xl)*(thetay2**x2)*(thetay3**x3) 
e  =  (thetaxl**yl)*(thetax2**y2)*(thetax3**y3) 
cl  =  d*e 

d  =  fctl*(thetayl**x2)*(thetay2**xl)*(thetay3**x3) 
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e  =  (thetaxl**y2)*(thetax2**yl)*(thetax3**y3) 

d2^fctl*(thetayl**x2)*(thetay2**x3)*(thetay3**xl) 
e  =  (thetaxl**y2)*(thetax2**y3)*(thetax3**yl) 

d3^fctl*(thetayl**x3)*(thetay2**x2)*(thetay3**xl) 
e  =  (thetaxl**y3)*(thetax2**y2)*(thetax3**yl) 

d  =fct  l6*  (thetay  1  **xl )  *  (thetay2**x3)  *  (thetay3**x2) 
e  =  (thetaxl**yl)*(thetax2**y3)*(thetax3**y2) 

d  ^fctie*(thetayl**x3)*(thetay2**xl)*(thetay3**x2) 

e  =  (thetaxl**y3)*(thetax2**yl)*(thetax3**y2) 
c6  =  d*e 

fct3  =  cl+c2+c3+c4+c5+c6 

prb  =  prb  +  fct3 

endif 

powr  =  powr  +  prb 

print*, j,powr 

end  do 

print* 

end 

subroutine  combo  (nn,rr, choose) 
double  precision  nn,rr,cc,diff,choose,dd,i 
diff  =  nn  -  rr 

if  (rr  .EQ.  0.0D0  .OR.  rr  .EQ.  nn)  then 

choose  =  1.0D0 
else 

if  (diff  .GE.  rr)  then 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,diff+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,rr 

dd  =  dd*i 
end  do 

choose  =  cc/dd 
else 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,rr-f  1.0DO,-1.0DO 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,diff 

dd  =  dd*i 
end  do 

choose  =  cc/dd 


endif 
endif 
return 
end 
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program  bnd2by2by2 

integer  i,j,k,s,mmm,totcnt,mm,m 

integer  c(100),cnt,critcnt,ii,jj 

double  precision  xll,x21,yll,y21,yyl(10000000) 

double  precision  tl,t2,xxl(10000000),hl,cl,c2 

double  precision  n,n2,xx2(10000000),c3,a,b 

double  precision  vl,v2,numer,mh,den,al,a2,c4,c5 

double  precision  yy2(10000000),c6,c7,xxlow(2),mxfnct 

double  precision  ntl,nt2,prbl,prb2,c8,bnd,delta(20) 

double  precision  pl,p2,p3,p4,ntl,nt2,xlow,ylow 

double  precision  xhigh,yhigh,prbl,prb2,pbndl,c62 

double  precision  pbnd2,pbnd3,pbnd4,bndprb,bndprb2 

double  precision  bound(100000),bound2(100000),mxbnd 

double  precision  lowl(8,10000),low2(8,10000),yylow(2) 

read*,n,bnd 

n2  =  2.0D0*n 

delta(l)  =  1.0D0/(2.0D0*n2) 

critcnt  =  0 
cnt  =  1 
j=0 
s  =  0 
k  =  0 

C 

C    Determine  critical  region 


c 


do  xll  =  0.0D0,n 
do  yll  =  0.0D0,xll 
do  x21  =  0.0D0,n 
do  y21  =  0.0D0,x21 

if  (xll  .EQ.  yll  .AND.  x21  .EQ.  y21)  then 

goto  10 

endif 

tl  =  xll  +  yll 

t2  =  x21  +  y21 

al  =  abs((xll-yll)/2.0D0) 

a2  =  abs((x21-y21)/2.0D0) 

numer  =  (al  +  a2  -  0.50D0)**2.0D0 

vl  =  (tl*(n2  -  tl))/(4.0D0*(n2-1.0D0)) 

v2  =  (t2*(n2  -  t2))/(4.0D0*(n2-1.0D0)) 

den  =  vl+v2 

mh  =  numer/den 

if  (mh  .GE.  bnd)  then 

critcnt  =  critcnt  +  1 

xxl  (critcnt)  =  xll 

xx2(critcnt)  =  x21 

yyl  (critcnt)  =  yll 

yy2(critcnt)  =  y21 


if  (xll  .NE.  yll  .AND.  x21  .NE.  y21)  then 
const  (critcnt)  =  4.0D0 

6ls6 

const  (critcnt)  =  2.0D0 

endif 
endif 
10  end  do 
end  do 
end  do 
end  do 

print*, "crit  done  !!" 
print*, critcnt 

C 

C    First  attempt  to  bound  function 
C 

mmm  =  0 
totcnt  =  0 
mm  =  0 
m  =  1 
c(l)  =  0 

hi  =  0.50D0  -  delta(m) 
print*, hi 

do  xlow  =  0.0D0,hl,delta(m) 
xhigh  =  xlow  +  delta(m) 
pi  =  xlow  +  delta(m)/2.0D0 
do  ylow  =  0.0D0,xlow,delta(m) 
yhigh  =  ylow  4-  delta(m) 
k  =  k  +1 

p2  =  ylow  +  delta(m)/2.0D0 
p3  =  1.0D0  -  pi 
p4  =  1.0D0  -  p2 
totcnt  =  totcnt  +  1 
bndprb  =  0.0D0 
bndprb2  =  0.0D0 
do  i  =  1, critcnt 
tl  =  xxl(i)+yyl(i) 
t2  =  xx2(i)+yy2(i) 

ntl  =  n2  -  tl 
nt2  =  n2  -  t2 
prbl  =  tl/n2 

prb2  =  t2/n2 

if  (prbl  .LE.  xlow)  then 

pbndl  =  xlow 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 


endif 
goto  20 

else 

pbndl  =  xhigh 
if  (prb2  .LE.  ylow)  then 
pbnd2  =  ylow 
else 

pbnd2  =  yhigh 

endif 
endif 

20  pbnd3  =  1.0D0  -  pbndl 
pbnd4  =  1.0D0  -  pbnd2 
call  combo(n,xxl(i),cl) 
call  combo(n,xx2(i),c2) 
call  combo(n,yyl(i),c3) 
call  combo(n,yy2(i),c4) 

c6  =  (pbndlC*3*tl)*(pbnd2**t2)*(pbnd3**ntl)*(pbnd4 
c62  =  (Pl**tl)*(p2**t2)*(p3**ntl)*(p4**nt2) 
c7  =  const(i)*c5*c6 
c8  =  const(i)*c5*c62 
bndprb  =  bndprb  +  c7 
bndprb2  =  bndprb2  +  c8 
end  do 

bndprb  =  bndprb 
bound2(k)  =  bndprb2 
if  (bndprb  .LE.  0.0250D0)  then 
mmm  =  mmm  +  1 
bound(mmm)  =  bndprb 
else 

c(cnt)  =  c(cnt)  +  1 

lowl(m,c(cnt))  =  xlow 

low2(m,c(cnt))  =  ylow 

endif 
end  do 
end  do 

print*,cnt,totcnt,mmm,c(cnt) 
if  (c(cnt)  .GT.  0)  then 
m  =  m  +  1 

delta(m)  =  delta(m-l)/2.0D0 

goto  30 

else 
goto  40 

endif 


* 


n 

C    Determining  bounds  for  the  above  regions  (  that  needed 
C    further  partitioning  )  iteratively,  until  the  desired 
C    accuracy  is  attained. 
C 

30  c(cnt+l)  =  0 
do  mm  =  l,c(cnt) 
xxlow(l)  =  lowl(m-l,mm) 
xxlow(2)  =  lowl(m-l,mm)  +  delta(m-l) 
yylow(l)  =  low2(m-l,mm) 
yylow(2)  =  low2(m-l,mm)  +  delta(m-l) 
do  ii  =  1,2 
do  jj  =  1,2 
xlow  =  xxlow(ii) 
xhigh  =  xlow  +  delta(m) 
pi  =  xlow  +  delta(m)/2.0D0 
ylow  =  yylow(jj) 
yhigh  =  ylow  +  delta(m) 
k  =  k  +1 

p2  =  ylow  +  delta(m)/2.0D0 

P3  =  1.0D0  -  pi 

p4  =  1.0D0  -  p2 

totcnt  =  totcnt  +  1 

bndprb  =  0.0D0 

bndprb2  =  0.0D0 

do  i  =  l,critcnt 

tl  =  xxl(i)+yyl(i) 

t2  =  xx2(i)+yy2(i) 

ntl  =  n2  -  tl 
nt2  =  n2  -  t2 
prbl  =  tl/n2 

prb2  =  t2/n2 

if  (prbl  .LE.  xlow)  then 

pbndl  =  xlow 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 

endif 
goto  20 

else 

pbndl  =  xhigh 

if  (prb2  .LE.  ylow)  then 

pbnd2  =  ylow 

else 

pbnd2  =  yhigh 
endif 
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25  pbnd3  =  l.ODO  -  pbndl 
pbnd4  =  l.ODO  -  pbnd2 
call  combo(n,xxl(i),cl) 
call  combo(n,xx2(i),c2) 
call  combo(n,yyl(i),c3) 
call  combo(n,yy2(i),c4) 
c5  _  ci*c2*c3*c4 

c6  =  (Pbndl**tl)*(pbnd2**t2)*(pbnd3**ntl)*(pbnd4**nt2) 

c62  =  (Pl**tl)*(p2**t2)*(p3**ntl)*(p4**nt2) 

c7  =  const(i)*c5*c6 

c8  =  const(i)*c5*c62 

bndprb  =  bndprb  +  c7 

bndprb2  =  bndprb2  +  c8 

end  do 

bndprb  =  bndprb 

bound2(k)  =  bndprb2 

if  (bndprb  .LE.  0.0250D0)  then 

mmm  =  mmm  +  1 

bound(mmm)  =  bndprb 

else 

c(cnt+l)  =  c(cnt+l)  +  1 

lowl(m,c(cnt+l))  =  xlow 

low2(m,c(cnt+l))  =  ylow 

endif 
end  do 
end  do 
end  do 

if  (c(cnt)  .GT.  0)  then 
m  =  m  +  1 

delta(m)  =  delta(m-l)/2.0D0 
print*,  m,delta(m) 
goto  30 
endif 

C 

C    Determining  accuracy  of  the  bound. 
C 

40  a  =  bound(l) 
do  j  =  2,mmm 
if  (bound(j)  .GT.  a)  then 
b  =  bound  (j) 

endif 

a  =  b 
end  do 

mxbnd  =  a 
a  =  bound2(l) 
do  j  =  2,k 


if  (bound2(j)  .GT.  a)  then 
b  =  bound2(j) 

endif 

a  =  b 
end  do 

mxfnct  —  a 

print*,mxbnd,mxbnd-mxfnct 
end 

subroutine  combo  (nn,rr,choose) 
double  precision  nn,rr,cc,difF,choose,dd, 
diff  =  nn  -  rr 

if  (rr  .EQ.  0.0D0  .OR.  rr  .EQ.  nn)  then 

choose  =  1.0D0 
else 

if  (diff  .GE.  rr)  then 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,diff+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,rr 

dd  =  dd*i 
end  do 

choose  =  cc/dd 
else 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,rr+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,diff 

dd  =  dd*i 
end  do 

choose  =  cc/dd 

endif 
endif 
return 
end 


program  fet2by2by2crit 
character  *9  fill,fil2 
character  *3  mh,numbr 
ch3.r3.ctcr  ^4  extn 

integer  i,j,k,s,cr,ncr,ns,l,m,ml,last,m2,kk,ll,ntemp 

double  precision  xll,x21,yll,y21,sxll(100000),nsx21(100000) 

double  precision  sx21(  100000)  ,inc,sp(  1000), tot,tl,t2 

double  precision  n,nl,ul,u2,perm,prob,probstar 

double  precision  dl,d2,d,nml,nm2,nm3,nm4,nm,n2,ttl 

double  precision  syll (100000), sy21  (100000), tt2,p,p2 

double  precision  sprob(  100000), nsprob(  100000), filler 

double  precision  nsyll(100000),nsy21(100000),nsxll(100000) 

alpha  =  0.05D0 

do  kk  =  2,4 

do  11  =  0,9 

mh  =  'cmh' 

numbr  =  'nmb' 

extn  =  '.dat' 

fill  =  mh  //  char(48+kk)  //  char (48+11)  //  extn 
fil2  =  numbr  //  char(48+kk)  //  char (48+11)  //  extn 
print*, fill 

ntemp  =  (10*kk)  +  11 
n  =  dble(ntemp) 
print*, n 
m  =  0 
1  =  0 

n2  =  2.0D0*n 
inc  =  0.80D0 
sp(l)  =  0.0D0 

last  =  idnint(n2)  +  1 
sp(last)  =  1000.0D0 
do  1  =  2,last-l 
sp(l)  =  sp(l-l)  +  inc 

end  do 

do  ml  =  2, last 

do  m2  =  2, last 

i  =  0 
j  =  0 


ns  =  0 

ncr  =  0 

tot  =  0.0D0 

nl  =  int(n/2) 

do  xll  =  0,nl 

do  yll  =  xll,n-xll 

do  x21  =  0,nl 


do  y21  =  x21,n-x21 

i  =  i  +  1      tl  =  xll  +  yll 

t2  =  x21  +  y21 

if  (tl  XT.  sp(ml)  .AND.  tl  .GE.  sp(ml-l))  then 

filler  =  0.02D0 

else 

goto  10 

endif  ,       v.  , 

if  (t2  .LT.  sp(m2)  .AND.  t2  .GE.  sp(m2-l))  then 

filler  =  0.02D0 

else 

goto  10 

endif 

call  combo(n2,tl,dl) 
call  combo(n2,t2,d2) 
d  =  dl*d2 

call  combo(n,xll,nml) 
call  combo(n,yll,nm2) 
call  combo(n,x21,nm3) 
call  combo(n,y21,nm4) 
nm  =  nml*nm2*nm3*nm4 
prob  =  nm/d 

ul  =  1.0D0 
u2  =  1.0D0 

if  (xll  .NE.  yll)  then 

ul  =  2.0D0 
endif 

if  (x21  .NE.  y21)  then 

u2  =  2.0D0 
endif 

perm  =  ul*u2 
probstar  =  perm*  prob 

if  (probstar  .LE.  alpha  .AND.  tl  .EQ.  t2)  then 

s  =  s  +  1 

sxll(s)  =  xll 

syll(s)  =  yll 

sx21(s)  =  x21 

sy21(s)  =  y21 

sprob(s)  =  probstar 

endif 

if  (probstar  .LE.  alpha  .AND.  tl  .NE.  t2)  then 

ns  =  ns  +  1 

nsxll(ns)  =  xll 

nsyll(ns)  =  yll 

nsx21(ns)  =  x21 

nsy21(ns)  =  y21 

nsprob(ns)  —  probstar 


endif 
10  end  do 
end  do 
end  do 

end  do      do  j  =  l,s 

pvalue  =  0.0D0 

p  =  sprob(j) 

tl  =  sxll(j)  +  syll(j) 

t2  =  8x210)  +  sy21(j) 

do  k  =  l,s 

ttl  =  sxll(k)  +  syll(k) 
tt2  =  sx21  (k)  +  sy21(k) 
p2  =  sprob(k) 

if  (tl  .EQ.  ttl  .AND.  t2  .EQ.  tt2)  then 

if  (p2  .LE.  p)  then 

pvalue  =  pvalue  +  p2 

end.if 
endif 
end  do 

if  (pvalue  .LE.  alpha)  then 
m  =  m  +  1 
open(l,FILE=fill) 

write(l,'(4fl2.5)')  sxll(j),syll(j),sx21(j),sy21(j) 

endif 
end  do 

do  j  =  l,ns 

pvalue  =  0.0D0 

p  =  nsprob(j) 

tl  =  nsxll(j)  +  nsyll(j) 

t2  =  nsx21(j)  +  nsy21(j) 

do  k  =  l,ns 

ttl  =  nsxll(k)  +  nsyll(k) 
tt2  =  nsx21(k)  +  nsy21(k) 
p2  =  nsprob(k) 

if  (tl  .EQ.  ttl  .AND.  t2  .EQ.  tt2)  then 

if  (p2  .LE.  p)  then 

pvalue  =  pvalue  +  p2 

end  if 

endif 

end  do 

if  (pvalue  .LE.  alpha)  then 
m  =  m  +  1 
open(l,FILE=fill) 

write(l,'(4fl2.5)')  nsxll(j),nsyll(j),nsx21(j),nsy21(j) 

endif 
end  do 
end  do 


end  do 

print*, m 

open(l,FILE=fil2) 

write(l,'(i25)')  m 

end  do 
end  do 

end  subroutine  combo  (nn,rr,choose) 
double  precision  nn,rr,cc,diff,choose,dd,i 
diff  =  nn  -  rr 

if  (rr  .EQ.  0.0D0  .OR.  rr  .EQ.  nn)  then 

choose  =  1.0D0 
else 

if  (diff  .GE.  rr)  then 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,diff+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,rr 

dd  =  dd*i 
end  do 

choose  =  cc/dd 
else 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,rr+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  1.0D0,diff 

dd  =  dd*i 
end  do 

choose  =  cc/dd 

endif 
endif 
return 
end 


program  fet2by2by2pwr 
integer  m,c,ul,u2,perm,m,tot 

double  precision  xl,x2,yl,y2,thetaxl,thetax2,thetayl,thetay2 

double  precision  n,nxl,nyl,cxl,cx2,cyl,cy2,fctl,al,a2,tabll5 

double  precision  tabl9,tabll0,tablll,tabll2,tabll3,tabll4,tabll6 

double  precision  tabll,tabl2,tabl3,tabl4,tabl5,tabl6,tabl7,tabl8 

double  precision  nx2,ny2,prb,powr,nthetaxl,nthetax2 

double  precision  nthetayl,nthetay2,prbl,prb2 

read*,n,thetaxl,thetax2,thetayl,thetay2 

powr  =  0.0D0 

c  =  0 
m  =  0 

open(l,FILE  =  'nmb.dat') 
read(l,'(i25)')  tot 
print*, tot 

open(l,FILE  =  'pwr.dat') 
do  m  =  l,tot 
c  =  c  +  1 

read(l,'(4fl2.5)')  xl,yl,x2,y2 

ul  =  1 
u2  =  1 
tl  =  xl  +  yl 
t2  =  x2  +  y2 
nxl  =  n  -  xl 
nyl  =  n  -  yl 
nx2  =  n  -  x2 
ny2  =  n  -  y2 

if  (xl  .NE.  yl  .AND.  tl  .NE.  n)  then 

ul  =  4 
endif 

if  (xl  .NE.  yl  .AND.  tl  .EQ.  n)  then 

ul  =  2 
endif 

if  (xl  .EQ.  yl  .AND.  xl  .NE.  nxl)  then 

ul  =  2 
endif 

if  (x2  .NE.  y2  .AND.  t2  .NE.  n)  then 

u2  =  4 

endif 

if  (x2  .NE.  y2  .AND.  t2  .EQ.  n)  then 

u2  =  2 
endif 

if  (x2  .EQ.  y2  .AND.  x2  .NE.  nx2)  then 

u2  =  2 

endif 

perm  —  ul*u2 

call  combo(n,xl,cxl) 

call  combo(n,x2,cx2) 


call  combo(n,yl,cyl) 
call  combo(n,y2,cy2) 

nthetaxl  =  l.ODO  -  thetaxl 
nthetax2  =  l.ODO  -  thetax2 
nthetayl  =  l.ODO  -  thetayl 
nthetay2  =  l.ODO  -  thetay2 
fctl  =  cxl*cx2*cyl*cy2 
if  (perm  .EQ.  16)  then 
al  =  (thetaxl**xl)*(nthetaxl**nxl^ 
a2  =  (thetax2**x2)*(nthetax2**nx2 

tabll  =  fctl*al*a2 
al  =  (thetaxl**xl)*(nthetaxl**nxl 

a2  =  (thetax2**y2)*(nthetax2**ny2 

tabl2  -  fctl*al*a2 
al  =  (thetaxl**yl)*(nthetaxl**nyl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tabl3  =  fctl*al*a2 
al  =  (thetaxl**yl)*(nthetaxl**nyl 

a2  =  (thetax2**y2)*(nthetax2**ny2 

tabl4  =  fctl*al*a2 
al  =  (thetaxl**xl)*(nthetaxl**nxl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabl5  =  fctl*al*a2 
al  =  (thetaxl**xl)*(nthetaxl**nxl 

a2  =  (thetax2**ny2)*(nthetax2**y2 

tabl6  =  fctl*al*a2 
al  =  (thetaxl**yl)*(nthetaxl**nyl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabl7  =  fctl*al*a2 
al  =  (thetaxl**yl)*(nthetaxl**nyl 

a2  =  (thetax2**ny2)*(nthetax2**y2 

tabl8  =  fctl*al*a2 
al  =  (thetaxl**nxl)*(nthetaxl**xl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tabl9  =  fctl*al*a2 
al  =  (thetaxl**nxl)*(nthetaxl**xl 

a2  =  (thetax2**y2)*(nthetax2**ny2 

tabllO  =  fctl*al*a2 
al  =  (thetaxl**nyl)*(nthetaxl**yl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tablll  =  fctl*al*a2 
al  =  (thetaxl**nyl)*(nthetaxl**yl 

a2  =  (thetax2**y2)*(nthetax2**ny2 

tabll2  =  fctl*al*a2 
al  =  (thetaxl**nxl)*(nthetaxl**xl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabll  3  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)* 


thetayl**yl)*(nthetayl**nyl 
thetay2**y2)*(nthetay2**ny2 

thetayl**yl)*(nthetayl**nyl 
thetay2**x2)*(nthetay2**nx2 

thetayl  **xl)*  (nthetayl  **nxl 
thetay2**y2)*(nthetay2**ny2 

thetayl**xl)*(nthetayl**nxl 
thetay2**x2)*(nthetay2**nx2 

thetayl**yl)*(nthetayl**nyl 
thetay2**ny2)*(nthetay2**y2 

thetayl**yl)*(nthetayl**nyl 
thetay2**nx2)*(nthetay2**x2 

thetayl**xl)*(nthetayl**nxl 
thetay2**ny2)*(nthetay2**y2 

thetayl**xl)*(nthetayl**nxl 
thetay2**nx2)*(nthetay2**x2 

thetayl**nyl)*(nthetayl**yl 
thetay2**y2)*(nthetay2**ny2 

thetayl**nyl)*(nthetayl**yl 
thetay2**x2)*(nthetay2**nx2 

thetayl**nxl)*(nthetayl**xl 
thetay2**y2)*(nthetay2**ny2 

thetayl**nxl)*(nthetayl**xl 
thetay2**x2)*(nthetay2**nx2 

thetayl**nyl)*(nthetayl**yl 
thetay2**ny2)*(nthetay2**y2 


(thetayl**nyl)*(nthetayl**yl) 


c 


a2  =  (thetax2**ny2)*(nthetax2**y2)*(thetay2**nx2)*(nthetay2**x2) 

alb=  tthetSl^ny 1)*  (nthetaxl  **y  1 )  *  (thetay  1  **nxl)*(nthetayl**xl) 
a2  =  (thetax2**nx2)*(nthetax2**x2)*(thetay2**ny2)*(nthetay2**y2) 
tabll5  =  fctl*al*a2 

al  =  (thetaxl**nyl)*(nthetaxl**yl)*(thetayl**nxl)*(nthetayl**xl) 

a2  =  (thetax2**ny2)  *  (nthetax2**y2)  *  (thetay2**nx2)  *  (nthetay2**x2) 

tabll6  =  fctl*al*a2  ,  tiA  ,  ir  ..0 
prbl  =  tabll+tabl2+tabl3+tabl4+tabl5+tabl6+tabl7+tabl8 

prb2  =  tabl9+tabll0+tablll+tabll2+tabll3+tabll4+tabll5+tabll6 

prb  =  prbl+prb2 

endif 

if  (perm  .EQ.  8  .AND.  xl  .EQ.  nyl)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl4  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**nx2)*(nthetax2**x2)*(thetay2**ny2)*(nthetay2**y2) 
tabl5  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**ny2)*(nthetax2**y2)*(thetay2**nx2)*(nthetay2**x2) 
tabl6  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**nx2)  *  (nthetax2*  *x2)  *  (thetay2*  *ny2)  *  (nthetay2**y2) 
tabl7  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**ny2)*(nthetax2**y2)*(thetay2**nx2)*(nthetay2**x2) 
tabl8  =  fctl*al*a2 

prb  =  tabll+tabl2+tabl3+tabl4+tabl5+tabl6+tabl7+tabl8 
endif 

if  (perm  .EQ.  8  .AND.  xl  .EQ.  yl)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)  *  (nthetax2**nx2)*  (thetay2**y2)  *  (nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)  *(nthetax2**ny2)  *  (thetay2**x2)  *  (nthetay2**nx2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 


c 


a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl4  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**nx2)*(nthetax2**x2)*(thetay2**ny2)*(nthetay2**y2) 
tabl5  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2*  *ny2)  *  (nthetax2*  *y2)  *  (thetay2*  *nx2)  *  (nthetay2**x2) 

tabl6  =  fctl*al*a2  ^, 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 

a2  =  (thetax2**nx2)  *  (nthetax2**x2)  *  (thetay2**ny2)  *  (nthetay2**y2) 

tabl7  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**ny2)*(nthetax2**y2)*(thetay2**nx2)*(nthetay2**x2) 
tabl8  =  fctl*al*a2 

prb  =  tabll+tabl2+tabl3+tabl4+tabl5+tabl6+tabl7+tabl8 
endif 

if  (perm  .EQ.  8  .AND.  x2  .EQ.  ny2)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2*  *x2)  *  (nthetax2**nx2)  *  (thetay2**y2)  *  (nthetay2*  *ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)*  (nthetax2**ny2)  *  (thetay2**x2)  *  (nthetay2**nx2) 
tabl'2  =  fctl*al*a,;) 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl4  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabl5  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**y2)  *  (nthetax2**ny2)  *  (thetay2**x2)*  (nthetay2**nx2) 
tabl6  =  fctl*al*a2 

al  =  (thetaxl**nyl)*(nthetaxl**yl)*(thetayl**nxl)*(nthetayl**xl) 
a2  =  (thetax2**x2)  *  (nthetax2**nx2)  *  (thetay2**y2)  *  (nthetay2**ny2) 
tabl7  =  fctl*al*a2 

al  =  (thetaxl**nyl)*(nthetaxl**yl)*(thetayl**nxl)*(nthetayl**xl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl8  =  fctl*al*a2 

prb  =  tabll+tabl2+tabl3+tabl4+tabl5+tabl6+tabl7+tabl8 
endif 
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if  (perm  .EQ.  8  .AND.  x2  .EQ.  y2)  then 


C 
C 


C 


thetayl**yl)*(nthetayl**nyl) 
thetay2**y2)*(nthetay2**ny2) 


al  =  (thetaxl**xl)*(nthetaxl**nxl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tabll  =  fctl*al*a2 
al  =  (thetaxl**xl)*(nthetaxl**nxl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabl2  =  fctl*al*a2 
al  =  (thetaxl**yl)*(nthetaxl**nyl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tabl3  =  fctl*al*a2 
al  =  (thetaxl**yl)*(nthetaxl**nyl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabl4  =  fctl*al*a2 
al  =  (thetaxl**nxl)*(nthetaxl**xl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tabl5  =  fctl*al*a2 
al  =  (thetaxl**nxl)*(nthetaxl**xl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabl6  =  fctl*al*a2 
al  =  (thetaxl**nyl)*(nthetaxl**yl 

a2  =  (thetax2**x2)*(nthetax2**nx2 

tabl7  =  fctl*al*a2 
al  =  (thetaxl**nyl)*(nthetaxl**yl 

a2  =  (thetax2**nx2)*(nthetax2**x2 

tabl8  =  fctl*al*a2 

prb  =  tabll +tabl2+tabl3+tabl4+tabl5+tabl6+tabl7+tabl8 
endif 


thetayl**yl)*(nthetayl**nyl) 
thetay2**ny2)*(nthetay2**y2) 

thetayl**xl)*(nthetayl**nxl) 
thetay2**y2)  *  (nthetay2*  *ny2) 

thetayl**xl)*(nthetayl**nxl) 
thetay2**ny2)*(nthetay2**y2) 

thetayl**nyl)*(nthetayl**yl) 
thetay2**y2)*(nthetay2**ny2) 

thetayl**nyl)*(nthetayl**yl) 
thetay2**ny2)*(nthetay2**y2) 

thetayl**nxl)*(nthetayl**xl) 
thetay2**y2)  *  (nthetay2**ny2) 

thetayl**nxl)*(nthetayl**xl) 
thetay2**ny2)*(nthetay2**y2) 


if  (ul  .EQ.  1  .AND.  u2  .EQ.  4)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)  *  (nthetax2**nx2)  *  (thetay2**y2)  *  (nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**nx2)*(nthetax2**x2)*(thetay2**ny2)*(nthetay2**y2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)  *  (nthetax2**ny2)*  (thetay2**x2)  *  (nthetay2**nx2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**ny2)*(nthetax2**y2)*(thetay2**nx2)*(nthetay2**x2) 

tabl4  =  fctl*al*a2 

prb  =  tabll  +  tabl2  +  tabl3  +  tabl4 

endif 

if  (ul  .EQ.  4  .AND.  u2  .EQ.  1)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 


a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**vl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**x2)  *  (nthetax2**nx2)  *  (thetay2**y2)  *  (nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**nyl)*(nthetaxl**yl)*(thetayl**nxl)*(nthetayl**xl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 

tabl4  =  fctl*al*a2 

prb  =  tabll  +  tabl2  +  tabl3  +  tabl4 

endif 

if  (ul  .EQ.  2  .AND.  u2  .EQ.2)  then 
if  (xl  .NE.  yl  .AND.  x2  .NE.  y2)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**y2)*  (nthetax2**ny2)  *  (thetay2**x2)  *  (nthetay2**nx2) 

tabl4  =  fctl*al*a2 

prb  =  tabll  +  tabl2  +  tabl3  +  tabl4 

endif 

if  (xl  .EQ.  yl  .AND.  x2  .NE.  y2)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**x2)  *  (nthetax2**nx2)  *  (thetay2**y2)  *(nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**nxl)*(nthetaxl**xl)*(thetayl**nyl)*(nthetayl**yl) 
a2  =  (thetax2**y2)  *  (nthetax2**ny2)  *  (thetay2**x2)  *  (nthetay2*  *nx2) 

tabl4  =  fctl*al*a2 

prb  =  tabll  +  tabl2  +  tabl3  +  tabl4 

endif 

if  (xl  .NE.  yl  .AND.  x2  .EQ.  y2)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 


^1 


tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**nx2)*(nthetax2**x2)*(thetay2**ny2)*(nthetay2**y2) 
tabl2  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**x2)  *  (nthetax2*  *nx2)  *  (thetay2*  *y2)  *  (nthetay2**ny2) 
tabl3  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**nx2)*(nthetax2**x2)*(thetay2**ny2)*(nthetay2**y2) 

tabl4  =  fctl*al*a2 

prb  =  tabll  +  tabl2  +  tabl3  +  tabl4 

endif 

endif 

if  (perm  .EQ.  2  .AND.  xl  .EQ.  yl)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2*  *x2)  *  (nthetax2**nx2)  *  (thetay2**y2)  *  (nthetay2*  *ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**y2)*(nthetax2**ny2)*(thetay2**x2)*(nthetay2**nx2) 

tabl2  =  fctl*al*a2 
prb  =  tabll  +  tabl2 
endif 

if  (perm  .EQ.  2  .AND.  x2  .EQ.  y2)  then 

al  =  (thetaxl**xl)*(nthetaxl**nxl)*(thetayl**yl)*(nthetayl**nyl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 
tabll  =  fctl*al*a2 

al  =  (thetaxl**yl)*(nthetaxl**nyl)*(thetayl**xl)*(nthetayl**nxl) 
a2  =  (thetax2**x2)*(nthetax2**nx2)*(thetay2**y2)*(nthetay2**ny2) 

tabl2  =  fctl*al*a2 
prb  =  tabll  +  tabl2 
endif 

powr  =  powr  +  prb 

end  do 

print* 

print*, powr 

end 

subroutine  combo  (nn,rr, choose) 

double  precision  nn,rr,cc,diff,choose,dd,i 

diff  =  nn  -  rr 

if  (rr  .EQ.  0.0D0  .OR.  rr  .EQ.  nn)  then 

choose  =  1.0D0 
else 

if  (diff  .GE.  rr)  then 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,diff+1.0D0,-1.0D0 
cc  ~  cc*i 


end  do 

do  i  =  1.0D0,rr 

dd  =  dd*i 
end  do 

choose  =  cc/dd 
else 

dd  =  1.0D0 
cc  =  1.0D0 

do  i  =  nn,rr+1.0D0,-1.0D0 

cc  =  cc*i 
end  do 

do  i  =  l.ODO.diff 

dd  =  dd*i 
end  do 

choose  =  cc/dd 

endif 
endif 

return 
end 
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